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EXECUTIVE  SUMMARY 


A  waterway  simulation  model  has  been  developed  to 
estimate  tow  delays  at  series  of  locks.  This  model  can  also 
estimate  tow  travel  time  along  waterways,  the  means  and 
variances  of  interarrival  and  interdeparture  times  at  each 
lock,  and  the  coal  inventory  levels  and  expected  stock-out 
amounts.  This  is  a  stochastic,  microscopic  and  event - 
scanning  model.  Each  simulation  run  takes  from  several 
seconds  to  several  minutes  on  a  PS/2  personal  computer 
depending  on  system  complexity. 

This  model  can  handle  any  distributions  for  trip 
generation,  lock  service  times  and  tow  sizes.  These 
distributions  can  be  specified  with  tables  to  best  represent 
reality  or  with  equations  to  represent  standard  statistical 
distributions.  Currently,  travel  speeds  are  assumed  to  be 
normally  distributed  while  general  distributions  based  on 
empirical  observations  are  used  for  trip  generation,  lock 
service  times  and  tow  sizes.  The  service  discipline  in  this 
model  is  FIFO.  The  model  ce s  accommodate  parallel  servers 
with  unequal  service  rates.  In  addition,  it  can  evaluate 
stall  effects  explicitly. 

Validation  results  show  that  the  model  works 
satisfactorily.  To  check  the  logic  of  this  model,  it  is 
compared  against  the  well  established  M/G/1  queues  over  a 
wide  range  of  volume  to  capacity  ratios  (0.04  to  0.89) .  The 
average  absolute  error  is  0.54%.  This  simulation  results 
are  also  compared  with  the  observed  data  at  five  lock  sites 
on  the  Mississippi  River.  The  validation  criteria  include 
average  waiting  times,  chamber  volumes,  and  cut  volumes. 

The  validation  results  show  that  the  simulation  model 
represents  the  real  system  quite  well. 

Although  this  simulation  model  requires  only  a  few 
seconds  to  a  few  minutes  for  each  lock  and  each  run  on  a 
PS/2  personal  computer,  that  is  still  hardly  affordable  for 
direct  application  in  large  combinatorial  network  investment 
problems . 

A  numerical  method  has  been  developed  for  estimating 
delays  through  a  series  of  G/G/1  queues  with  inflows  and 
outflows  occurring  only  at  end  nodes.  This  numerical  method 
is  an  approximation  of  the  simulation  model.  It  can  quickly 
evaluate  large  combinatorial  investment  problems  and  thus 
helps  identify  the  best  combinations  of  investment 
alternatives  for  further  analysis  (e.g.,  by  simulation). 

This  numerical  method  was  originally  developed  for 
systems  with  bi-directional  servers.  With  a  few 
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simplifications,  this  method  can  be  adapted  for  the  m^ re 
generally  env^-^antered  systems  with  one -directional  servers. 
The  two-way  algorithm  employs  an  iterative  alternating 
direction  scanning  procedure  to  estimate  the  interarrival 
and  interdeparture  time  distributions  lock  by  lock  until  the 
interdeparture  time  variances  for  successive  iterations 
converge.  The  performance  of  this  two-way  algorithm  is 
tested  with  satisfactory  results.  The  one-way  algorithm 
only  scans  the  interarrival  time  and  interc eparture  time 
distributions  from  the  firs'"  to  the  last  lock  without  any 
iteration  and  should,  theoretically,  be  less  subject  to 
interdependence  errors. 

Both  the  two-way  and  one-way  algorithms  rely  on  several 
metamodels  estimated  from  the  results  of  a  previously 
developed  simulation  model  These  metamodels  provide  the 
following  valuable  results  for  series  of  G/G/j.  queues. 

1.  The  delay  metamodel  (Eq.  41)  indicates  how  the  V/C 
ratios,  interarrival  time  distributions, 
inteidep -■rture  time  distributions  and  service 
distributions  affect  the  average  waiting  times  for 
G/G/1  queues.  This  delay  function  is  an  exact 
solution  based  on  Marshall's  formula  for  the 
variance  of  interdeparture  times. 

2.  The  relations  among  the  coefficients  of  variation 
of  interdeparture  times,  interarrival  tires, 
service  times,  and  the  V/C  ratio  are  formulated  in 
the  departures  metamodel .  The  structure  of  the 
departure  function  (Eq.  36)  is  based  on  functions 
for  the  squared  coefficients  of  variation  of 
interdeparture  times.  By  applying  Laplace 
transforms,  these  functions  (i.e.,  Eqs .  31  and  32) 
are  derived  theoretically  in  this  study. 
Statistical  estimation  of  the  parameters  yields  a 
very  good  fit.  The  function's  standard  error  of 
0.0058  is  extremely  tight  compared  to  its  mean  of 
0.8311.  The  parameters  also  have  very  tight 
standard  errors.  In  addition,  this  departure 
function  is  consistent  with  Burke's  Theorem.  The 
results  show  that  the  metamodeling  approach 
combining  queuing  theory  and  statistical 
estimation  based  on  simulation  outputs  is  quite 
successful.  This  approach  for  approximating 
departure  processes  snould  be  very  useful  for 
analyzing  networks  of  queues. 

3.  The  arrivals  module  provides  the  relation  between 
the  variance  of  interarrival  times  and  the 
variance  of  interdeparture  times  from  the  adjacent 
queue  stations  when  speed  variations  change  the 
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headway  distributions  between  successive  queue 
stations . 

It  should  be  noted  that  the  metamodel  parameters  have 
been  estimated  from  data  bases  representing  conditiop=  on 
waterways.  However,  results  to  date  indicate  that  such  good 
metamodels  can  be  extended  to  other  applications. 

The  applications  of  the  numerical  method  are  currently 
limited  to  series  of  G/G/1  queues  in  waterways  with  inflows 
and  outflows  occurring  at  end  nodes.  The  following 
extensions  are  desirable  to  increase  its  applicability. 

1.  A  function  should  be  developed  to  estimate  the 
variance  of  interarrival  times  when  inflows  and 
outflows  occur  between  queuing  stations.  Such  an 
arrival  function  should  reflect  the  complex 
superposition  and  splitting  of  flows  in  general 
networks  of  queues,  including  tree  and  grid 
networks.  A  proposed  approach  would  compute  the 
variance  of  overall  arrival  rates  as  the  sum  of 
the  arrival  rate  variances  from  all  inflows, 
assuming  individual  inflows  are  independent  of 
each  other,  and  then  develop  the  relation  between 
the  variance  of  interarrival  times  and  the 
variance  of  arrival  rates. 

2.  The  effects  of  random  failures  (i.e.,  stalls) 
might  be  incorporated  by  treating  stalls  as  a 
second  class  of  users  with  its  own  arrival  and 
service  time  distributions.  The  present  method 
can  already  incorporate  stalls  as  part  of  the 
exogenously  specified  service  time  distribution, 
to  the  extent  that  stalls  are  related  to  traffic 
volumes.  However,  it  is  desirable  to  have  a 
method  which  can  evaluate  stall  effects 
explicitly . 

3.  If  possible,  the  numerical  method  should  be 
extended  to  locks  with  two  or  more  dissimilar 
chambers.  This  is  rather  difficult  for  chambers 
with  different  characteristics,  because  the 
chamber  assignment  process  affects  lock  capacity. 

4 .  The  ability  to  handle  unequal  directional  trip 
rates  would  be  desirable  for  many  applications, 
although  waterways  usually  operate  with  roughly 
equal  directional  flows. 

5.  It  would  be  desirable  to  model  queues  with  limited 
storage  space,  which  are  highly  unusual  on 
waterways  but  fairly  common  in  dense  road 
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networks,  computers,  and  other  queuing  network 
applications . 

Additional  statistical  and  computational  tests  are  also 
desirable  to  further  validate  this  method  and  to  extend  its 
applicability . 

The  approach  developed  in  the  numerical  method  may  be 
applied  not  only  to  a  lock  queuing  system,  but  also  to  some 
other  systems  of  queues.  More  comprehensive  simulation 
experiments  are  desirable  for  developing  the  arrival  and 
departure  process  modules  for  networks  of  G/G/k  queues. 

These  simulation  experiments  may  include  wider  ranges  of 
variables,  and  other  special  distributions  for  interarrival 
and  service  time  distributions.  Further  research  in 
estimating  interarrival  time  distributions  with  multiple 
unequal  inflows  would  be  necessary  for  extending  the 
numerical  method  to  general  networks  of  queues  which  may 
have  inflows  and  outflows  at  any  node. 

The  final  methodology  for  estimating  waterway  delays 
may  combine  simulation  and  numerical  method.  The  numerical 
method  may  be  used  in  analyzing  large  combinatorial 
problems,  that  may  be  encountered  in  investment  scheduling 
or  real  time  traffic  control.  The  numerical  method  may  be 
used  to  screen  alternatives  so  that  only  a  few  of  the  most 
promising  ones  are  examined  more  thoroughly  with  the 
simulation  model.  Guidelines  should  be  developed  for 
switching  between  the  numerical  method  and  simulation  in 
applications  requiring  intermediate  accuracy. 

The  feasibility  of  approximating  simulation  models  for 
complex  queuing  systems  with  simple  metamodels  should  be 
intensively  explored.  If  feasible,  it  would  have  important 
applications  in  combinatorial  optimization  and  real-time 
control  for  transportation  and  communication  networks, 
computers,  and  manufacturing  systems. 
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CHAPTER  1 
INTRODUCTION 


Backgrovind 

Inland  waterway  transportation  is  quite  important  in 
the  U.S.  and  other  countries,  especially  for  heavy  or  bulky 
commodities,  since  it  is  inexpensive,  energy  efficient,  and 
safe.  Most  U.S.  waterways  consist  of  stepped  navigable 
pools  formed  by  dams  across  natural  rivers.  The  lock 
structures  used  to  raise  or  lower  vessels  between  adjacent 
pools  constitute  the  major  bottlenecks  in  the  U.S.  waterway 
network  [27]  and  generate  extensive  queues.  The  resulting 
delays  and  variability  of  service  times  have  very 
substantial  economic  implications. 

Some  locks  have  only  one  chamber,  while  others  may  have 
two  parallel  chambers  whose  characteristics  may  differ.  The 
most  common  chamber  sizes  are  110'  x  1200'  (i.e.,  110  feet 

wide  and  1200  feet  long)  and  110'  x  600'.  Each  chamber  size 
can  accommodate  a  limited  number  of  barges  at  one  time.  For 
example,  a  110'  x  1200'  chamber  can  accommodate  at  most  17 
standard  barges  plus  a  towboat,  while  a  110'  x  600'  chamber 
can  accommodate  at  most  8  standard  barges  plus  a  towboat . 

If  a  tow  has  more  barges  than  the  chamber  can  accommodate, 
it  must  be  disassembled  into  several  pieces  (called  "cuts") 
to  move  through  the  chamber  and  must  later  be  reassembled. 
Therefore,  the  service  time  distributions  depend  on  chamber- 
size  and  tow-size  distributions.  Sometimes,  chambers  will 
be  out  of  service  (i.e.,  "stalled")  for  various  reasons  such 
as  freezing,  accidents,  and  mechanical  failures. 

Figure  1  shows  a  simple  diagram  of  a  lock  queuing 
system.  Locks  are  the  servers  and  tows  are  customers 
waiting  to  be  served  by  locks.  In  the  lock  queuing  system, 
tows  from  both  directions,  upstream  and  downstream,  share 
the  same  lock  servers,  while  in  most  other  queuing  systems 
the  servers  are  exclusively  one-directional.  In  this  work, 
the  term  "two-way  traffic  operations"  characterizes  the  lock 
queuing  system  while  "one-way  traffic  operations" 
characterizes  queuing  systems  in  which  servers  are 
exclusively  one-directional. 

Arrival-time  and  service-time  distributions  at  locks 
are  fairly  complex.  Carroll  [5]  and  Desai  [14]  found  that 
service  times  are  not  exponentially  distributed,  and 
arrivals  are  not  Poisson  distributed.  Other  standard 
distributions  have  been  tested  for  the  present  study, 
without  consistent  success.  Locks  with  a  single  chamber  may 
be  modeled  as  G/G/l  queuing  systems.  (The  notation  means 
"generally  distributed  arrivals/  generally  distributed 
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FIGURE  1,  LOCK  QUEUING  SYSTEM 


service  times/  one  server".)  However,  locks  with  two 
parallel  chambers  may  not  be  treated  simply  as  G/G/2  queuing 
systems  unless  the  parallel  chambers  are  identical. 

The  lock  service-time  distributions  are  affected  by  the 
chamber  assignment  discipline  at  locks  with  two  dissimilar 
chambers.  There  the  "main"  chamber  is  larger  than  the 
"auxiliary"  chamber  and  can  accommodate  without  disassembly 
large  tows  that  might  require  several  cuts  and  far  larger 
service  times  to  move  through  the  auxiliary  chamber. 

However,  if  the  same  number  of  cuts  is  required  through 
either  chamber,  a  shorter  auxiliary  chamber  may  provide 
faster  service  since  entry  and  exit  times  may  be  reduced. 
Therefore,  lock  service-time  distributions  are  dynamic  and 
depend  on  the  chamber  assignment  discipline. 

Considerable  interdependence  may  exist  among  locks  in  a 
series.  The  departure  distributions  differ  from  the  arrival 
distributions  since  the  service-time  distributions  change 
the  tow  headways.  The  departures  from  one  lock  usually 
affect  the  arrivals  at  the  next  lock.  Therefore,  it  is 
improper  to  assume  that  the  locks  are  independent .  The 
interdependence  among  locks  increases  the  difficulty  in 
estimating  delays  for  the  lock  queuing  system  since  it  is 
necessary  at  each  lock  to  identify  the  interarrival -time 
distributions  of  flows  from  adjacent  locks. 

Two-way  traffic  operation  through  common  servers 
complicates  the  interdependence  of  lock  delays  and  precludes 
the  use  of  some  otherwise  interesting  queuing  models. 

Delays  are  determined  by  the  arrival  distributions  and 
service-time  distributions.  It  is  much  more  difficult  to 
identify  the  arrival  distributions  for  two-way  traffic 
systems  than  for  one-way  traffic  systems.  The  arrival 
distribution  at  one  lock  is  affected  by  departures  from  both 
upstream  and  downstream  locks,  while  departures  from  this 
lock  also  affect  the  arrivals  at  upstream  and  downstream 
locks.  ^or  example,  in  Figure  1  the  arrivals  at  Lock  2 
would  be  affected  by  the  departures  from  Locks  1  and  3 .  The 
departures  from  Locks  1  and  3  toward  Lock  2  are  highly 
correlated  with  the  arrivals  at  1  and  3  from  2 .  Some  of  the 
arrivals  at  1  and  3  represent  departures  from  2.  Hence,  the 
arrival  distributions  of  these  three  locks  are 
interdependent.  Thus,  two-way  traffic  operation  complicates 
the  estimation  of  the  two  arrival  distributions  at  each 
lock.  The  arrival  distributions  depend  not  only  on  the 
departure  distributions  from  the  adjacent  locks,  but  also  on 
speed  variations  and  distances  between  adjacent  locks. 

Random  failures,  which  in  inland  waterways  are  called 
stalls,  contribute  significantly  to  the  difficulties  in 
estimating  delays.  Stalls,  which  interrupt  lock  operations 
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and  thereby  increase  delays  and  service  time  variances,  are 
relatively  rare  compared  to  other  events,  however  their 
occurrence  is  very  difficult  to  predict. 

Problem  Statement 

A  reliable  and  efficient  method  for  estimating  delays 
is  essential  for  evaluating  and  scheduling  waterway 
investments.  When  many  investment  proposals  are  considered, 
their  selection  and  scheduling  becomes  a  large  combinatorial 
problem,  since  waterway  investments  are  often 
interdependent . 

The  purpose  of  this  research  is  to  develop  an  efficient 
and  reliable  method  to  estimate  delays  through  a  series  of 
waterway  locks.  The  difficulties  in  estimating  delays  for 
such  a  lock  queuing  system  are  summarized  as  follows: 

1.  Arrival-time  and  service-time  distributions  are 
generally  distributed  (i.e.,  they  cannot  be 
represented  by  any  standard  statistical 
distribution) . 

2.  The  service  rates  of  parallel  chambers  may  be  very 
different . 

3.  Service-time  distributions  are  affected  by  the 
chamber  assignment  discipline. 

4  Considerable  interdependence  exists  among  locks  in 
series . 

5.  Two-way  traffic  operation  through  bi-directional 
chambers  complicates  the  analysis. 

6 .  The  arrival  distributions  depend  not  only  on  the 
departures  from  previous  locks  but  also  on  the 
distances  and  speed  distributions  between  locks. 

7.  Stalls  increase  tne  means  and  variances  of  delays. 

Even  neglecting  the  special  complexities  associated 
with  waterways,  the  available  analytic  sclutions  for 
estimating  delays  through  a  series  of  G/G/l  queues  are  quite 
inadequate.  These  analytic  solutions  are  approximations  or 
are  subject  to  some  limitations.  They  will  be  discussed  in 
greater  detail  in  the  literature  review. 

Scope 


This  research  seeks  to  estimate  delays  along  a  series 
of  locks  in  waterways.  The  arriva^.s  ^nd  service  times  at 


4 


these  locks  are  generally  distributed.  Some  locks  may  have 
two  parallel  unequal  chambers.  All  chambers  must  handle 
upstream  and  downstream  traffic.  Also,  these  locks  are 
subject  to  interruptions  of  operation  due  to  stalls. 

This  research  also  seeks  to  explore  the  feasibility  of 
substituting  single  equations  for  simulations  of  complex 
systems  of  queues.  If  feasible,  such  substitutes  for 
simulation  would  have  important  applications  in 
combinatorial  optimization  and  real-time  control  for  various 
transportation  and  communication  networks,  computers,  and 
manufacturing  systems. 
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CHAPTER  2 
LITERATURE  REVIEW 


The  literature  reviewed  in  this  study  concentrates  on 
lock  delay  models,  lock  interdependence,  waterway  simulation 
models,  single  queues,  and  networks  of  queues  which  are  the 
subjects  of  greatest  relevance  to  this  work.  Of  these,  the 
work  on  waterway  simulation  models  and  on  networks  of  queues 
represent  the  closest  parallels  to  the  models  developed  for 
this  study. 

Lock  Delay  Models 

Two  models  based  on  the  application  of  queuing  theory 
have  been  found  for  estimating  lock  delays.  DeSalvo  and 
Lave  [13]  represent  the  lock  operation  as  a  simple  single¬ 
server  queuing  process  with  Poisson  distributed  arrivals  and 
exponentially  distributed  service  times.  Wilson  [37] 
extends  the  previous  model  by  treating  the  service  processes 
as  general  distributions.  Both  models  are  designed  for 
analyzing  single-lock  delays. 

Lave  and  DeSalvo  [13]  proposed  that  approximate  lock 
delays  could  be  estimated  with  an  M/M/1  queuing  model 
(Poisson  arrivals/Exponential  service  times/1  server) . 
Unfortunately,  their  assumptions  about  the  arrival  and 
service-time  distributions  do  not  satisfactorily  fit  the 
physical  system  of  locks  cn  waterways.  In  particular, 
Carroll  and  Wilson  [7]  found  that  the  exponentially 
distributed  service  times  do  not  correspond  well  with 
empirical  evidence.  Therefore,  it  is  not  proper  to  assume 
that  the  service  times  are  exponentially  distributed. 

The  assumption  of  Poisson  distributed  arrivals  does  not 
fit  every  lock  in  the  waterway  system.  Carroll  and  Desai 
[5,14]  studied  the  arrival  processes  for  40  locks  on  the 
Illinois,  Mississippi,  and  Ohio  River  systems,  utilizing 
1968  data.  The  results  of  Chi-square  tests  show  that  13  out 
of  40  locks  had  non-Poisson  arrivals  at  the  5%  significance 
level.  Therefore,  the  assumption  of  Poisson  distributed 
arrivals  is  questionable. 

Wilson  [37]  proposed  an  M/G/l  model  (Poisson 
arrivals/general  service  times/1  server)  to  analyze  lock 
delays.  Unlike  DeSalvo  and  Lave  [13],  Wilson  represented 
the  service  processes  as  generally  distributed  rather  than 
exponentially  distributed,  which  is  far  more  realistic  [7] . 
However,  even  the  Poisson  arrivals  assumption  is  not 
realistic  for  all  locks.  Another  deficiency  in  Wilson's 
model  is  that  no  exact  queuing  results  are  available  for 
locks  with  two  chambers  in  parallel.  Thus,  Wilson's  model 
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can  only  be  applied  to  locks  with  Poisson  distributed 
arrivals  and  a  single  chamber. 

Two  other  deficiencies  exist  in  both  of  the  above 
models.  First,  neither  of  these  models  accounts  for  stalls. 
Stalls  cause  service  interruptions  at  locks,  thus  reducing 
lock  capacities  or  increasing  delays.  Their  occurrence  is 
very  difficult  to  predict.  Thus,  Kelejian's  efforts  to 
model  stalls  and  stall  durations  have  not  yet  yielded  strong 
results  despite  the  rigorous  statistical  methods  employed 
[19] .  The  second  deficiency  is  that  both  models  were 
developed  to  analyze  delays  at  a  single  lock.  Since  the 
delays  at  adjacent  locks  may  be  highly  related,  it  is 
desirable  to  analyze  lock  delays  for  entire  systems. 

Based  on  the  above  discussion,  it  seems  desirable  to 
develop  a  model  that  can  represent  the  physical  conditions 
satisfactorily  and  analyze  systems  of  locks.  The  physical 
conditions  that  should  be  represented  include  generally 
distributed  arrivals  and  service  times,  and  service 
interruptions  due  to  stalls. 

Lock  Interdependence 

The  departure  process  of  a  queue  in  a  network  is  of 
special  interest  because  it  is  likely  to  determine  the 
arrival  process  at  the  following  queue  in  the  network.  This 
produces  interdependence  among  system  elements.  In  other 
words,  lock  interdependence  can  be  expressed  in  terms  of  the 
relations  among  the  arrival,  service,  and  departure 
processes  at  one  lock.  If  the  distributions  of  tow  arrivals 
and  departures  at  one  lock  are  identical,  it  can  be  shown 
[11]  that  this  lock  is  independent  of  other  locks. 

Burke's  theorem  [11]  states  that  the  steady-state 
output  of  a  stable  M/M/m  queue  with  input  parameter  X  and 
service  time  parameter  fx  for  each  of  the  m  channels  is  in 
fact  a  Poisson  process  with  the  same  rate  X.  Therefore,  for 
two  locks  in  series,  if  the  arrival  process  is  Poisson 
distributed  and  the  service  process  is  exponentially 
distributed  at  the  first  lock,  the  second  lock  will  be 
independent  of  the  first  lock.  In  addition,  the  second  lock 
will  have  the  same  arrival  distribution  as  the  first  lock. 

Burke's  theorem  established  that  an  M/M/m  queue  is 
independent  of  the  other  processes  in  a  network  of  queues . 
Hence,  the  decomposition  technique  can  be  applied  to  analyze 
the  delays  of  this  network  (if  all  queues  in  it  are  M/M/m) . 
While  this  is  a  powerful  result,  the  delays  at  inland 
waterway  locks  cannot  be  analyzed  independently  since  the 
service  processes  are  observed  to  be  non-exponentially 
distributed  [7] . 
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Carroll  and  Desai  [5,14]  studied  tow  arrival  processes 
at  locks  on  the  inland  waterway  system.  They  assumed  that, 
if  the  tow  arrivals  at  a  lock  follow  a  Poisson  distribution, 
the  queuing  processes  at  this  lock  are  independent  of  what 
has  occurred  at  the  preceding  lock.  They  employed  a 
waterway- system  simulation  model  to  analyze  operations 
through  a  series  of  locks.  The  goodness-of - f it  of  the  tow 
interarrival  distributions  to  a  negative  exponential 
distribution  was  then  tested  at  all  locks,  using  the  Chi- 
Square  test.  The  results  showed  that  13  of  40  locks  had 
non-Poisson  arrivals  at  the  5%  significance  level.  Hence, 
Poisson  distributed  arrivals  may  not  be  assumed  for  the 
inland  waterway  system. 

Waterway  Simulation  Models 

The  system  simulation  models  developed  to  analyze  lock 
delays  and  tow  travel  times  originated  mainly  from  Howe's 
microscopic  model  [17] .  In  that  model  service  times  were 
based  on  em.pirically-determined  frequency  distributions.  To 
avoid  some  troublesome  problems  and  errors  associated  with 
the  requirement  to  balance  long-run  flows  in  Howe's  model, 
Carroll  and  Bronzini  [6]  developed  another  waterway- system 
simulation  model.  It  provided  detailed  outputs  on  such 
variables  as  tow  traffic  volumes,  delays,  processing  times, 
transit  times,  averages  and  standard  deviations  of  delay  and 
transit  times,  queue  lengths,  and  lock  utilization  ratios. 

Both  of  the  above  models  simulate  waterway  operations 
in  detail  but  require  considerable  amounts  of  data  and 
computer  time,  which  limit  their  applicability  for  problems 
with  large  networks  and  numerous  combinations  of  improvement 
alternatives.  Both  models  assume  Poisson  distributions  for 
tow- trip  generation,  which  is  not  always  realistic.  More 
importantly  for  reliability  analyses,  neither  of  these 
models  explicitly  accounts  for  stalls,  which  are  very 
different  in  frequency  and  duration  from  other  events  and 
have  significant  effects  on  overall  transit-time 
reliability . 

Hence  a  waterway  simulation  model  that  explicitly 
accounts  for  stalls  is  desirable  for  evaluating  and 
scheduling  lock  improvement  projects. 

Single  Queues 

Queuing  theory  deals  with  the  stochastic  behavior  of 
service  systems.  Historically,  queuing  theory  was  first 
developed  in  the  context  of  telephone  traffic.  Erlang,  in 
particular,  made  many  important  contributions  to  the  subject 
in  the  early  part  of  this  century.  Telephone  traffic 
remained  one  of  the  principal  applications  until  about  1950. 
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Recently,  queuing  theory  has  been  applied  to  many  other 
fields  including  communications,  computer  networks, 
transportation,  and  manufacturing. 


Elementary  queuing  theory  is  based  on  pure  Markov 
queues.  In  a  Markov  process,  the  state  of  the  queuing 
system  at  any  time  depends  only  on  the  current  state  and  not 
on  any  previous  states.  A  very  import" ant  special  class  of 
Markov  processes  is  the  birth-death  process  in  which  state 
transitions  take  place  between  neighboring  states  only.  The 
application  of  the  birth-death  process  has  helped  in 
developing  exact  equilibrium  solutions  for  queues  with 
Poisson  distributed  arrival  processes  and  exponentially 
distributed  service  times. 


G/M/1  or  M/G/l  queues  are  non-Markovian  stochastic 
processes  that  m.ay  be  analyzed  with  one  of  the  following 
four  approaches.  First,  the  method  of  imbedded  Markov 
chains  focuses  on  the  number  of  customers  present  in  a 
queuing  system  immediately  following  a  departure  [21] . 
Second,  the  method  of  stages,  proposed  by  Erlang,  has  the 
disadvantage  that  it  merely  gives  a  procedure  to  find  the 
solution  but  does  not  show  the  solution  as  an  explicit 
expression  [21] .  Third,  G/G/1  systems  are  solved  using 
Lindley's  integral  equation.  Therefore,  this  approach  can 
be  applied  to  the  special  cases  of  G/M/1  or  M/G/i  [21] . 
Fourth,  the  method  of  supplementary  variables  uses  the  state 
vector  [N(t),Xg(t)]  where  N(t)  is  the  number  of  customers  in 
the  system  at  time  t  and  XgCt)  is  the  service  time  already 
experienced  by  the  customer  served  at  time  t  [21] . 


G/G/1  queues  are  difficult  cases  in  queuing  theory  and 
the  available  techniques  for  handling  them  arc  incomplete. 

As  mentioned  above,  Lindley's  integral  eq”ation  is 
applicable  for  G/G/1  systems.  The  waiting-time  distribution 
for  customers  W(y)  in  the  system  G/G/1  can  be  written  as 


W  (y)  = 


where 


fy 

W(y-x) dC(x) 

J  -00 
1 


if  y  2  0 
if  y  <  0 


(1) 


y  =  waiting  time 

X  =  the  difference  between  the  service  time  for  one 
customer  and  the  interarrival  time  between  this 
customer  and  the  next . 

C(x)  =  the  probability  distribution  function  (PDF)  of  u 


In 


addition, 
w  =  ^ 


the  mean  waiting  time 

+  +  E [ta] ^  (1  -  p) ^ 

2E  [t;^]  (1-p) 


is 


E  [I^] 
2E[I] 


(2) 


10 


where 


Op^-.  variance  of  the  interarrival  times 
a^-.  variance  of  service  times 
t;^:  interarrival  time 
p:  volume  to  capacity  ratio 
I  =  idle  period 

Unfortunately,  it  is  very  difficult  to  solve  the  two  moments 
(E[I^]  and  E[I])  of  the  idle  period  since  this  period 
depends  on  the  particular  way  in  which  the  previous  busy 
period  terminated. 

Solving  G/G/m  queues  is  even  more  difficult  than 
solving  G/G/1  queues.  The  methods  of  approximations  and 
bounds  have  been  proposed  to  solve  G/G/m  queues  [28,29] . 
These  are  accurate  and  efficient  under  heavy- traffic 
conditions . 


Bertsimas  [4]  proposed  a  methodology  for  analyzing  the 
waiting  times  of  G/G/m  queues  with  First-In-First-Out  (FIFO) 
service  discipline  and  mixed  generalized  Erlang  distributed 
arrivals  and  service  times.  This  method  could  be  applied  to 
more  realistic  situations  than  Poisson  arrivals  and 
exponential  service  times.  However,  without  a  departure 
function,  this  result  is  difficult  to  extend  to  a  series  of 
locks . 


When  traffic  is  heavy,  and  thus  the  ratio  p  approaches 
1,  the  average  waiting  time  for  G/G/1  queues  can  be  written 
as 


W 


2(1  -  p)E[t;^] 


(3) 


There  are  two  types  of  approximation  methods  for  G/G/m 
queues  [28,29]  .  First,  the  fluid  approximation  uses  the 
mean  values  to  represent  the  stochastic  processes.  Second, 
the  diffusion  approximation  improves  the  fluid  approximation 
by  describing  the  processes  with  means  given  by  the  fluid 
approximation  but  with  a  normal  distribution  describing  the 
fluctuations  about  that  mean.  The  important  assumption 
behind  these  two  approximation  methods  is  that  the  queue 
never  empties,  which  is  consistent  with  the  heavy  traffic 
assumption . 


Wolff  [38]  compared  the  delays  under  different  orders 
of  service  (service  disciplines)  and  concluded  that  the 
order  of  service  would  not  change  the  mean  of  waiting  times, 
but  the  service  discipline  Last - In- First -Out  (LIFO)  has  a 
larger  variance  than  the  First- In-First -Out  (p’IFO) 
discipline . 
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Networks  of  Queues 


Modeling  networks  of  queues  is  very  useful  for 
analyzing  the  performance  of  complex  systems,  such  as 
communication,  computer,  transportation,  and  manufacturing 
networks.  Exact  solutions  are  only  available  for  the 
Markovian  networks  of  queues  [2] .  For  more  realistic 
networks  of  queues  whose  arrivals  are  not  Poisson  and  whose 
service  times  are  not  exponential,  approximation  methods  are 
usually  employed  foi  system  performance  analysis. 

In  a  Markovian  network  of  queues,  all  arrival  processes 
in  the  network  are  Poisson  distributed  and  all  service  times 
are  exponentiallv  distributed.  The  queue-length 
distribution  of  individual  queues  in  the  Markovian  networks 
is  the  same  as  that  of  a  queue  with  Poisson  arrival  and 
departure  processes  [2] .  In  addition,  the  joint  queue 
length  distribution  is  the  product  of  the  individual  queue- 
length  distributions  [3,8,9,16,18,23].  That  is,  the  joint 
queue-length  distribution  has  the  product  form.  Therefore, 
the  decomposition  technique  can  be  applied  to  the  Markovian 
networks  of  queues  and  each  individual  queue  is  analyzed 
independently . 

The  concept  of  decomposition  can  also  be  applied  to 
more  general  (non-Markovian  and  one-way)  networks  of  queues. 
The  general  networks  of  queues  are  decomposed  approximately 
into  individual  queues.  These  individual  queues  are 
analyzed  independently  and  then  the  results  are  recombined 
for  evaluating  an  entire  system.  The  Queuing  Network 
Analyzer  (QNA)  [35]  ,  based  on  decomposition,  is  the  one  mo.st 
relevant  to  this  study.  The  decomposition  technique  is 
applied  in  other  examples  [10,22,30,31]. 

QNA  is  a  comprehensive  software  package  for 
approximating  the  congestion  measures  of  open  networks  of 
queues  with  multiple  servers.  The  service  discipline  in  QNA 
is  FIFO  and  the  network  capacity  is  unlimited.  The  most 
important  feature  of  QNA  is  that  the  arrivals  need  not  be 
Poisson  distributed  and  the  service  times  need  not  be 
exponentially  distributed.  QNA  employs  two  parameters, 
namely  the  mean  and  the  coefficient  of  variation,  to 
approximately  characterize  the  arrival  processes  and  then 
analyzes  the  individual  queuing  nodes  independently.  The 
decomposition  procedures  in  QNA  include  three  sc  ages  of 
analyses:  superposition  of  arrival  streams,  departure 
processes,  and  splitting  of  departure  streams.  The 
superposition  of  arrival  processes  in  QNA  is  a  variation  of 
the  hybrid  method  in  Albin  [1] .  The  departure  processes  in 
QNA  are  identified  by  the  stationary- interval  method  [35] . 
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To  approximate  the  superposition  of  arrival  processes, 
Albin  [1]  and  Whitt  [34]  suggested  the  use  of  the  hybrid 
method,  which  is  a  combination  of  the  stationary- interval 
method  and  asymptotic  method.  The  hybrid  method  is  mainly 
used  for  estimating  the  coefficient  of  variation.  The 
hybrid  coefficient  of  variation  of  intervals  between  two 
consecutive  arrivals,  is  a  convex  combination  of  the 

squared  coefficients  of  variation  obtained  from  the 
stationary-interval  method  and  asymptotic  method: 

Ch^  =  +  (l-w)C2^  (4) 

where 

:  squared  coefficient  of  variation  obtained  from  the 
asymptotic  method  [34] 

C2^ :  squared  coefficient  of  variation  obtained  from  the 
stationary- interval  method  [34] 
w:  weighting  function,  0  s  w  s  1,  w  =  w(p,  n*) 

p:  volume  to  capacity  ratio  of  the  queue 

n* :  effective  number  of  component  processes 

n"* :  [(X]^^  +  ...  +  X^^)/X^]'^ 

X^:  the  rate  of  the  ith  component  process 

X :  X  =  Xj^  +  ...  +  Xj^ 

Albin  [1]  proposed  that  the  weighting  function  should 
approach  1  as  p  approaches  1  and  approach  0  as  n*  approaches 
infinity.  Based  on  such  a  conjecture,  Albin  created  a  list 
of  candidate  weighting  functions  and  used  simulation  to 
identify  the  best  one.  Thus,  the  weighting  function  was 
developed  empirically.  She  also  developed  a  different 
weighting  function  for  each  operating  characteristic  to  be 
approximated.  For  example,  she  suggested  using  w= [1+6(1- 
p) ^'^n*] '^for  approximating  the  expected  number  in  the 
system;  w= [1  +  8 . 3 ( 1 -p ) ^ ■  ^n*]  for  approximating  the  standard 
deviation  of  the  number  in  the  system;  and  w=0  for 
approximating  the  probability  of  delay. 

It  should  be  noted  that  Albin' s  study  [1]  mainly 
estimated  the  coefficient  of  variation  for  the  superposed 
arrival  processes  to  queues.  The  weighting  function  in  her 
study  is  a  function  of  the  congestion  level  (p)  and  varies 
for  evaluating  different  operating  characteristics  (the 
expected  number  in  the  system,  the  standard  deviation  of  the 
number  in  system,  and  the  probability  of  delay) .  However, 
theoretically,  for  a  single  queuing  station,  the  coefficient 
of  variation  of  the  superposed  arrival  processes  should  be 
independent  of  the  operating  characteristics  of  queues. 

Thus,  the  coefficient  of  variation  should  not  be  affected  by 
p  and  should  not  vary  for  different  evaluation  purposes. 
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To  understand  why  p  is  influential  in  Albin's  weighting 
function  we  need  to  know  how  she  obtained  the  results.  In 
fact,  the  weighting  function  in  Albin's  study  was  not 
developed  directly  based  on  the  coefficient  of  variation  but 
based  on  the  approximating  operating  characteristics  [1] . 

She  obtained  a  weighting  function  that  can  produce  close 
agreement  with  the  approximating  operating  characteristics. 
Therefore,  her  results  may  be  affected  by  the  approximating 
operating  characteristics  and  not  represent  the  true 
coefficient  of  variation.  Since  the  superposition  arrival 
processes  in  QNA  is  a  variation  of  Albin's  hybrid  method, 

QNA  also  has  the  deficiency  described  above. 

Departure  processes  are  very  important  for  networks  of 
queues  since  they  often  become  the  inputs  for  other  queues. 
Whitt  [36]  tried  four  methods,  the  asymptotic  method,  the 
stationary- interval  method,  the  lag-1  correlation,  and 
hybrid  method,  for  approximating  the  departure  processes  by 
characterizing  two  parameters  or  the  first  two  moments. 

The  asymptotic  method  tries  to  approximate  the 
parameters  of  departure  processes  by  matching  their  long-run 
behavior.  The  departure  process  approximated  by  the 
asymptotic  method  is  just  the  arrival  process.  Therefore, 
the  result  obtained  by  the  asymptotic  method  may  not  be 
useful  for  capturing  the  interdependence  between  queues. 

The  potential  drawback  of  the  stationary- interval 
method  is  that  it  does  not  take  account  of  the  dependence 
among  successive  intervals.  Whitt  [36],  based  on  Marshall's 
expression  for  variances  of  departure  processes  [26] , 
characterized  the  coefficient  of  variation  of  departures  in 
terms  of  the  waiting  time  as  follows: 

+  2p2cg2  _  2X(l-p)W  (5) 

where 

:  squared  coefficient  of  variation 
times 

:  squared  coefficient  of  variation 
times 

Cg^ :  squared  coefficient  of  variation 

p:  volume  to  capacity  ratio 

X :  average  arrival  rate 

W:  average  waiting  time 

However,  the  exact  average  waiting  time  for  G/G/1  queues  is 
not  available.  Kraemer  and  Langenbach-Belz  [35]  developed 
the  following  approximation  formula: 

W  =  p(Ca"  +  Cs2)g/2/x{l  -  p)  (6) 


of  interdeparture 
of  interarrival 
of  service  times 
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where 


/X :  the  average  service  rate 

g  =  g  (p  ,Cj^^  .Cs^) 

g(p,C^^,Cs^)  =  exp  [-2  (1-p)  (1-C;^^)  ^/3p  (C;^^  +  Cg^)],  <  1 

=  1,  ^  1 

Whitt  [36]  combines  the  above  approximation  while  g  is 
assumed  to  be  one.  The  expression  of  is  as  follows: 

=  p2cg2  +  (1  _  p2)C^2 

Therefore,  the  coefficient  of  variation  obtained  with  the 
stationary- interval  method  by  Whitt  is  based  on  two 
successive  approximations  and  may  generate  errors  when 
is  less  than  1,  which  is  very  common  in  waterway  queuing 
networks . 

The  lag-1  correlation  provides  a  compromise  between  the 
stationary- interval  method  and  the  asymptotic  method.  It 
tries  to  match  the  local  behavior  of  a  departure  process  and 
also  to  partially  account  for  the  dependence  among 
successive  intervals.  However,  Whitt's  results  did  not  show 
any  improvement  when  the  lag-1  correlation  method  was  used. 

The  hybrid  method  for  the  departure  process  is  also  a 
convex  combination  of  the  parameters  from  the  asymptotic 
method  and  the  stationary- interval  method.  However,  this 
approach  shows  no  improvement  either  compared  to  the 
stationary- interval  method  [36]  . 

Whitt  [36]  concluded  that  the  stationary- interval 
method  has  the  best  performance  among  these  four  methods. 
Therefore,  the  departure  process  in  QNA  uses  the  stationary- 
interval  approach.  However,  Whitt's  approach  [36] 
approximates  the  departure  process  twice  and  may  generate 
significant  error  when  is  less  than  one.  Therefore,  it 
is  desirable  to  develop  a  method  to  approximate  the 
departure  process  more  precisely. 

Albin  [2]  also  tried  to  approximate  the  departure 
process  for  a  single  server  with  exponentially  distributed 
service  times  by  using  a  hybrid  approach.  The  building 
blocks  in  her  hybrid  method  are  the  asymptotic  method  and 
the  Poisson  method.  The  Poisson  method  approximates  the 
departure  process  as  a  Poisson  process.  The  asymptotic 
method  approximates  the  departure  process  as  the  arrival 
process.  Therefore,  the  in  Albin  [2]  becomes: 

Cd^  =  +  (1  _  w)  (8) 
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where 


w:  weighting  function,  w  s  vi  {n*  ,  p-^ ,  P2) 

n* :  the  effective  number  of  component  arrival 

processes 

p^:  the  volume  to  capacity  ratio  at  queuing  station  i 

It  is  notable  that  Albin's  approach  for  approximating 
departure  processes  has  the  same  deficiency  associated  with 
superposition  arrival  processes  since  the  weighting  function 
depends  on  the  volume  to  capacity  ratio  at  the  second 
queuing  station.  Theoretically,  the  departure  processes  at 
the  first  queuing  station  should  not  be  affected  by  the 
operation  at  the  second  queuing  station  unless  there  is  a 
spillback  from  the  second  queuing  station  that  blocks  the 
operation  of  the  first  queuing  station.  Also,  the  departure 
processes  should  not  vary  for  different  evaluation  purposes. 
In  addition,  Albin's  results  are  only  applicable  to  queues 
with  exponential  service  times.  Therefore,  it  is  desirable 
to  develop  a  methodology  to  accurately  estimate  departure 
processes  for  generally  distributed  arrivals  and  service 
times . 
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CHAPTER  3 

MICROSCOPIC  SIMULATION  MODEL 


A  simulation  model  has  been  developed  to  analyze  tow 
operations  along  waterways.  It  may  be  used  to  determine  the 
relations  among  delays,  tow  trips,  distributions  of 
generated  tow  trips,  lock  operations,  lock  service-time 
distributions,  travel  times,  coal  consumption,  and  coal 
inventories.  The  simulation  model  can  take  into  account 
stochastic  effects  and  seasonal  variations.  This  model 
enables  the  following  to  be  estimated:  tow  delays  at  each 
lock,  interarrival  and  interdeparture- time  distributions  for 
each  lock  and  for  each  direction,  tow  travel  times  along  che 
waterway,  inventory  levels  and  expected  stock-out  amounts 
for  coal  de'',ivered  by  waterway,  and  many  other  variables  of 
interest  to  waterway  users  and  managers.  The  estimation  of 
tow  delays,  tow  travel  times,  inventory  levels,  and  expected 
stock-out  amounts  is  useful  for  estimating  economic  benefits 
of  lock  improvements.  Moreover,  the  interarrival  and 
interdeparture  time  distributions  and  the  delays  estimated 
with  this  model  should  be  useful  for  analyzing  series  of 
queues  with  generally-distributed  interarrival  and  service¬ 
time  distributions. 

Data  Base 

The  simulation  model  is  developed  on  the  basis  of  PMS 
(lock  Performance  Monitoring  System)  data  collected  since 
1975.  This  data  base  includes  very  detailed  information  on 
traffic  through  the  locks  as  well  as  physical  aspects  of 
lockage  [15] .  It  is  very  useful  for  understanding  and 
quantifying  waterway  characteristics,  such  as  lock 
operations,  arrival  distributions,  service-time 
distributions,  tow-size  distributions,  and  stalls. 

Features  of  Simulation  Model 

The  simulation  model  developed  in  the  early  stages  of 
the  study  focuses  on  the  relations  among  delays,  trip 
generations,  distributions  of  tow  speeds,  arrivals, 
departures  and  service  t.’mes,  and  coal  inventories.  The 
output  of  this  model  includes  delays,  means  and  variances  of 
interarrival- time  and  incerdeparture- time  distributions,  and 
inventory  level.  These  are  useful  for  estimating  inventory 
costs,  stock-out  costs,  and  expected  benefits  due  to  lock 
rehabilitation  or  lock  construction.  The  relations  among 
delays,  interarrival,  interdeparture  and  service- time 
distributions  provide  results  for  the  analysis  of  series  or 
networks  of  queues. 
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The  simulation  model  is  microscopic.  It  traces  the 
movement  of  each  individual  tow  and  records  its 
characteristics,  including  its  number  of  barges,  commodity 
types,  speed,  origin  and  destination,  travel  direction,  and 
arrival  time  at  various  points.  In  addition,  the  model 
determines  cumulative  deliveries,  cumulative  consumption, 
and  actual  inventories  at  various  plants. 

The  simulation  model  is  an  event -scanning  model  where 
the  system  status  is  updated  by  events.  The  system  status 
includes  the  simulation  clock,  lock  operating  condition, 
movement  of  tows,  and  inventory  level  at  each  utility  plant. 
The  heart  of  the  simulation  model  is  the  scheduler  routine. 
This  routine  provides  the  control  for  the  entire  length  of 
the  simulation  period.  The  scheduler,  at  appropriate  times, 
invokes  all  other  operational  routines  necessary  to  process 
the  simulation. 

This  model  can  handle  any  distributions  for  trip 
generation,  travel  speeds,  lock  service  times  and  tow  sizes. 
These  distributions  can  be  specified  for  each  interval  in 
tables  or  by  standard  statistical  distributions.  Currently, 
travel  speeds  are  assumed  to  be  normally  distributed,  while 
general  distributions  based  on  empirical  observations  are 
used  for  other  input  variables.  Tows  are  allowed  to 
overtake  other  tows.  A  FIFO  (First-In-First-Out)  service 
discipline  is  currently  employed.  This  model  simulates  two- 
way  traffic  through  common  servers  and  accounts  for  stalls. 

There  are  five  types  of  events  in  this  model.  First, 
tow  trips  are  generated  stochastically  based  on  the  actually 
observed  traffic  distributions.  Currently,  the  model  uses  a 
table  to  represent  the  trip-generation  pattern  and  is, 
therefore,  not  limited  to  standard  mathematical  probability 
distributions . 

Second,  the  tow  entrance  in  a  lock  is  determined  by  tow 
arrival  time  at  that  lock,  the  times  when  chambers  become 
available,  and  the  chamber  assignment  discipline.  If  a  tow 
arrives  before  the  lock  is  available,  it  needs  to  wait  in 
the  queue  storage  area.  Otherwise,  it  is  served  according 
to  the  chamber  assignment  discipline  that  will  be  discussed 
later.  In  general,  the  lock  service  is  presently  First-In- 
First  -Out,  subject  to  the  chamber  assignment  procedure. 

Third,  whenever  a  coal  tow  arrives  at  its  destination, 
the  cumulative  delivery  at  the  destination  will  be  increased 
by  the  amount  of  coal  that  tow  is  carrying.  The  cumulative 
consumption  and  inventory  at  the  destination  will  also  be 
updated  then. 


18 


Fourth,  the  status  of  cumulative  consumptions, 
inventories,  and  consumption  rates  for  all  coal  destinations 
will  be  updated  in  every  unit  of  time.  This  provides 
detailed  information  on  inventory  levels  for  all  coal 
destinations . 

Fifth,  whenever  a  stall  occurs  at  a  chamber,  the 
chamber  becomes  unavailable  until  the  end  of  the  stall. 

The  size  of  waterway  systems  that  the  model  can  handle 
is  limited  by  the  computer  capacity  and  the  storage  capacity 
of  the  Fortran  compiler  or  linker.  The  simulation  model  has 
been  developed  with  "dynamic  dimensioning"  to  the  degree 
allowed  by  the  computer  system  available.  Parameter 
statements  are  used  so  that  the  dimensions,  and  hence 
capacities,  of  the  model  components  may  readily  be  modified. 
This  allows  the  maximum  flexibility  of  waterway  system 
design  and  the  most  efficient  computer  utilization.  Thus, 
the  dynamic  dimensioning  programming  technique  allows 
flexibility  in  the  number  of  locks,  chambers,  cuts,  waterway 
links,  tows,  utility  plants,  origin-destination  (0-D)  pairs 
and  simulation  time  periods.  Currently  this  model  can 
simulate  two-way  operation  on  a  mainline  waterway. 

This  simulation  model  is  programmed  in  Fortran- 77, 
which  allows  us  to  simulate  relatively  complex  operations. 
The  following  is  a  more  detailed  description  of  how  tow  trip 
generation,  tow  travel  times,  and  coal  inventory  levels  are 
computed  in  this  model.  The  overall  structure  of  the 
simulation  model  is  displayed  in  Figure  2.  A  logic  flow 
chart  is  provided  in  Figure  3. 

Assumptions  of  Simulation  Model 

The  simulation  model  was  developed  based  on  the 
following  assumptions: 

1.  The  time  interval  between  successive  tow  trip 
generations,  service  times,  and  tow  sizes  are 
generally  distributed.  These  distributions  are 
represented  by  probability  distribution  tables 
which  are  exogenous  inputs.  Therefore,  the 
simulation  model  may  be  applied  to  any  systems  of 
lock  systems . 

2.  The  tow  maintains  a  constant  size  through  the 
entire  trip. 

3.  The  service  discipline  is  First - In- First -Out 
(FIFO) .  This  assumption  is  consistent  with 
waterway  operations  on  the  Mississippi  and  Ohio 
Rivers . 
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FIGURE  2.  STRUCTURE  AND  ELEMENTS  OF  SIMULATION  MODEL 
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FIGURE  3.  FLOW  CHART  OF  SIMULATION  MODEL 
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4 .  The  main  chamber  is  preferred  for  chamber 
assignment .  This  reflects  the  wish  of  lock 
operators  to  avoid  the  additional  work  and  delays 
due  to  disassembling  and  re -assembl ing  extra  cuts 
through  an  auxiliary  chamber. 

5.  The  queue  storage  area  is  unlimited.  Storage 
space  is  not  a  significant  problem  on  waterways. 
However,  the  assumption  requires  modification  when 
the  simulation  model  or  the  numerical  method 
(developed  based  on  the  results  of  the  simulation 
model)  is  applied  to  other  systems. 

6.  Tow  speeds  are  normally  distributed.  For  lac-; 
better  data,  this  assumption  is  based  on  wate.-  / 
statistics  of  vessel  performance  [32] . 

7.  Each  tow  maintains  a  constant  speed  between  its 
origin  and  destination. 

8.  The  average  backhaul  speed  is  a  constant  ratio  of 
the  average  linehaul  speed.  The  assumption 
reflects  the  effects  of  currents  and  barge  loads. 

9.  The  time  intervals  between  two  successive  stalls 
and  the  durations  of  stalls  are  exponentially 
distributed . 

10.  The  barge  payload  is  the  same  for  all  coal  barges 
and  remains  constant  throughout  a  trip.  This 
implies  that  the  industry  has  a  standard  barge 
size  and  loads  each  barge  to  capacity. 

11.  There  is  an  exogenously  specified  fraction  of 
barges  carrying  coal  on  a  coal  tow.  That  is,  a 
coal  tow  also  carries  other  commodities.  This 
assumption  is  also  consistent  with  the  barge 
industry  operations. 

12.  The  consumption  rates  at  each  utility  plant  are 
uniformly  distributed  within  specified  (and 
arbitrarily  short)  intervals. 

Simulation  Routines 

The  simulation  model  consists  of  five  operation 
routines  and  one  scheduler  routine.  The  operation  routines 
are  associated  with  the  five  types  of  events:  (1)  tow  trip 
generation  (trip  generation  routine),  (2)  tow  arrival  at  a 
lock  (lock  routine),  (3)  coal  tow  arrival  at  its  destination 
(node  routine),  (4)  periodical  inventory  and  consumption 
updates  (periodical  routine),  and  (5)  stall  occurrence 
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(stall  routine) .  These  operation  routines  are  invoked  by 
the  scheduler  routine. 

Trip  Generation  Routine 

Tow  trips  are  generated  stochastically,  but  the  mean  of 
their  generating  distribution  is  constant  for  each  origin- 
destination  (0-D)  pair  over  each  simulation  time  period. 

Each  0-D  pair  has  its  own  distribution  for  trip  generation 
which  may  be  represented  in  this  model  by  a  probability 
distribution  table  or  a  standard  statistical  distribution. 
The  probability  distribution  tables  represent  cumulative 
distribution  curves,  where  the  abscissas  represent 
cumulative  frequency  and  the  ordinates  represent  the  ratio 
of  the  tabulated  variable  to  its  mean.  The  trip  generation 
tables  can  be  easily  changed  since  they  are  specified 
explicitly  in  the  input  data. 

Whenever  a  tow  trip  is  generated,  its  size  (numbers  of 
barges  per  tow)  and  speed  are  also  generated.  This  model 
assumes  that  each  tow  will  maintain  its  size  and  speed 
throughout  its  entire  trip.  As  in  trip  generation,  the  long 
run  average  tow  size  is  constant  for  each  0-D  pair,  but  the 
size  of  each  individual  tow  is  generated  stochastically. 

The  distribution  of  tow  sizes  is  represented  by  a 
probability  distribution  table  and  each  0-D  pair  has  its  own 
distribution  table.  The  tow  size  distribution  table  is  an 
exogenous  input  data.  Therefore,  it  flexibly  accommodates 
any  type  of  tow  size  distribution. 

Tow  speeds  are  specified  as  an  input  to  the  model  in 
the  form  of  a  probability  distribution.  For  lack  of  better 
data,  the  distribution  of  speeds  is  currently  assumed  to  be 
normal  [32] .  This  model  currently  assumes  that  tows 
maintain  constant  speeds  between  origins  and  destinations 
and  that  backhaul  speeds  are  a  constant  ratio  of  linehaul 
speeds . 

To  avoid  generating  extreme  speed  values,  a  speed  range 
is  specified.  Currently,  if  the  generated  speeds  are  lower 
than  the  2.5  percentile  speed  or  zero,  or  higher  than  the 
97.5  percentile  speed,  these  speeds  must  be  regenerated. 

Tow  traffic  is  divided  into  coal  traffic  and  non-coal 
traffic.  Therefore,  for  the  same  0-D  trips,  there  may  be 
two  different  0-D  pairs  with  different  trip  rates  for  coal 
traffic  and  for  non-coal  traffic,  respectively.  It  is 
presently  assumed  that  only  a  specified  fraction  of  the 
barges  on  a  coal  tow  are  carrying  coal.  ’  In  addition,  the 
randomly  generated  tow  sizes  for  coal  traffic  are  restricted 
by  upper  and  lower  limits  to  preclude  extremely  small  or 
extremely  large  tows.  The  tow  size,  the  fraction  of  barges 
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carrying  coal,  and  the  barge  payload  determine  the  amount  of 
coal  delivered  by  a  coal  tow. 

After  the  characteristics  associated  with  the  newly 
generated  tow  trip  are  determined,  the  trip  generation 
routine  also  determines  the  next  lock  the  tow  must  pass 
through  and  the  tow's  arrival  time  at  this  lock.  The  tow's 
arrival  time  at  this  lock  is  determined  by  its  trip 
generation  time  and  its  travel  time  over  the  distance 
between  the  origin  and  this  lock.  Meanwhile,  the  trip 
generation  routine  also  determines  the  a',  rival  sequence  at 
this  lock  for  the  tow  and  updates  the  next  tow's  arrival 
time  at  this  lock.  The  information  about  next  tow  arrival 
time  at  this  lock,  the  tow  arrival  sequence  at  this  lock, 
and  the  associated  tow  characteristics  (size,  speed,  travel 
direction)  is  recorded  for  reference  by  other  routines. 

Lock  Routine 


The  lock  routine  is  the  most  complex  one  in  the 
simulation  model.  It  performs  the  following  tasks: 

1.  chamber  assignment; 

2.  determination  of  tow  entrance  time; 

3.  computation  of  tow  waiting  time; 

4.  determination  of  chamber  available  time  (tow 
departure  time)  for  next  tow; 

5.  computation  of  interarrival  time  between  two 
successive  tows; 

6.  computation  of  interdeparture  time  between  two 
successive  tows; 

7.  determination  of  next  stop  (a  lock  or  the 
destination  node)  for  the  tow; 

8.  computation  of  travel  time  between  this  lock  and 
the  next  stop;  and, 

9.  determination  of  arrival  sequence  of  the  tow  at 
the  next  stop. 

The  lock  service  discipline  is  currently  First-In- 
First-Out  in  the  lock  routine.  The  lock  routine  serves  the 
first  tow  based  on  the  arrival  sequence.  Such  an  assumption 
simplifies  the  simulation  model  and  is  consistent  with 
waterway  operations  on  the  Mississippi  and  Ohio  Rivers. 
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The  chamber  assignment  is  based  on  the  chamber 
available  time,  the  expected  chamber  service  time,  and  the 
lock  selection  bias  factor.  The  chamber  available  time  is 
either  the  departui'e  time  of  previous  tow  or  the  stall 
ending  time.  The  expected  chamber  service  time  depends  on 
the  number  of  cuts  required  for  the  tow  moving  through 
chamber.  The  lock  selection  bias  factor  reflects  in 
equivalent  time  units  the  advantage  of  utilizing  the  main 
chamber . 

The  number  of  cuts  is  determined  by  cut  s^zes  nd  tow 
sizes.  The  cut  size  is  exogenous.  Each  chamber  has  an 
upper  limit  on  cut  size  (number  of  barges  that  can  be 
handled  simultaneously) ,  which  determines  the  number  of  cuts 
required.  A  tow  may  be  disassembled  into  different  numbers 
of  cuts  at  different  lock  chambers.  Therefore,  the  lock 
routine  needs  to  determine  the  required  number  of  cuts  at 
each  chamber  for  each  tow. 

It  is  noted  that  the  expected  chamber  service  time  used 
in  chamber  assignment  is  the  average  value  of  service  tin.c: 
for  the  required  cut  numbers.  This  reflects  the  reality  of 
operation.  In  fact,  the  lock  operator  would  not  know  in 
advance  the  actual  service  time  for  the  tow.  The  operator 
could  only  roughly  estimate  the  average  service  time  for  the 
tow  based  on  its  size. 

If  a  lock  has  dissimilar  chambers  in  parallel,  (main 
and  auxiliary  chambers  are  usually  provided) ,  it  is 
currently  assumed  that  the  main  chamber  will  be  preferred, 
unless  the  additional  wait  time  it  requires  (compared  to  the 
auxiliary  chamber)  exceeds  a  specified  level.  This  lock 
selection  bias  factor  reflects  the  additional  work  and 
delays  required  to  break  tows  into  more  (and  smaller)  cuts, 
move  them  separately  through  an  auxiliary  chamber  and  then 
reassemble  them.  Bias  factors  have  been  estimated 
separately  for  various  locks  from  empirical  data.  They  can 
be  easily  modified  since  they  are  specified  explicitly  in 
the  input  data. 

Once  the  chamber  is  assigned,  the  tow  entrance  time  is 
also  determined.  The  tow  entrance  time  is  either  the  tow 
arrival  time  or  the  available  chamber  time,  whichever  is 
greater.  If  the  tow  arrives  before  the  chamber  is 
available,  the  tow  must  wait  in  the  queue  storage  area. 
Otherwise,  the  tow  enters  its  assigned  chamber  immediately. 
The  queue  storage  area  is  currently  assumed  to  have 
unlimited  capacity  and  to  be  adjacent  to  the  lock. 

For  most  applications  the  tow  waiting  time  is  the  most 
important  output  of  the  simulation  model.  It  is  defi.ned  as 
the  difference  between  a  tow's  arrival  tim.e  and  its  entry 
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time  into  the  lock.  Such  tow  waiting  times  may  occur  well 
before  traffic  levels  approach  lock  capacity  since  tow 
arrivals  and  lock  service  times  are  not  constant. 

The  available  chamber  time  for  the  next  tow  is  the  tow 
departure  time  which  depends  on  tne  entrance  time  and  the 
lock  service  time  of  the  current  tow  and  on  the  possible 
occurrence  of  a  stall.  In  general,  without  a  stall 
occurrence,  the  tow  departure  time  is  defined  as  the  sum  of 
the  tow  entrance  time  and  the  required  lock  service  time. 

If  a  stall  occurs  before  the  service  is  completed,  the  tow 
would  be  detained  until  the  end  of  the  stall.  Therefore, 
the  tow  departure  time  is  equal  to  the  sum  of  the  tow 
entrance  time,  the  required  lock  service  time,  and  stall 
duration  if  there  is  a  stall. 

Lock  service  times  may  be  generated  from  a  specified 
distribution  table  or  a  standard  statistical  distribution. 
The  distribution  table  can  directly  reflect  actually 
observed  service  times.  Therefore,  the  model  can  be  applied 
to  any  types  of  locks.  Lock  service  times  will  be  affected 
by  lock  improvements  which  are  represented  by  smaller 
average  lock  service  times  or  reduced  service  time 
variances.  The  average  lock  service  tj.mes  vary  for 
different  locks,  chambers,  and  numbers  of  cuts. 

Whenever  a  stall  occurs  before  lock  service  is 
completed,  the  lock  routine  must  update  the  stall  event  at 
this  lock.  That  includes  information  on  the  stall 
occurrence  time  and  duration. 

The  interarrival  time  between  two  successive  tows  is 
defined  as  the  time  intervrl  between  the  current  tow's 
arrival  time  and  the  arrival  time  of  a  previous  tow. 
Similarly,  tiie  interdepar  ure  time  :  a  defined  as  the  time 
interval  between  the  current  tow's  departure  time  and  the 
departure  time  of  a  previous  tow. 

The  lock  routine  also  determines  the  next  stop  for  the 
tow  and  the  tow  arrival  time  at  the  next  stop.  The  next 
stop  could  be  a  lock  or  the  destination  of  the  tow.  The  tow 
arrival  time  at  the  next  stop  is  determined  by  the  tow 
departure  time  and  tow  travel  time  over  the  distance  between 
this  lock  and  the  next  stop.  If  the  next  stop  is  a 
destination  and  the  •  is  not  a  coal  tow,  the  lock  routine 
will  let  the  tow  le^  _  the  system;  otherwise,  the  lock 
routine  will  determine  the  tow's  arrival  sequence  at  the 
next  stop  and  update  the  next  tow's  arrival  time  ac  this 
stop  in  the  meantime.  The  information  about  the  next  tow's 
arrival  time  at  this  step,  the  tow  arrival  sequence  at  this 
lock,  and  the  associated  tow  characteristics  (size,  speed, 
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and  travel  direction)  is  recorded  for  reference  by  other 
routines . 

Node  Routine 


Whenever  a  coal  tow  arrives  at  its  destination,  i.e., 
at  a  utility  plant,  the  node  routine  updates  the  inventory 
status  there.  The  node  routine  computes  the  cumulative 
deliveries,  cumulative  consumption,  inventory  level,  and 
stock-out  amounts  and  durations. 

Cumulative  deliveries  are  determined  from  initial 
inventory  levels,  delivery  amounts,  and  arrival  time  at  the 
destination.  The  initial  inventory  level  is  an  exogenous 
parameter  specified  for  each  destination.  The  node  routine 
can  estimate  several  different  cumulative  deliveries 
associated  with  different  initial  inventory  levels 
simultaneously.  The  delivered  amount  is  determined  from  the 
barge  payload  and  the  number  of  arriving  coal  barges. 
Currently,  the  barge  payload  is  assumed  to  be  constant  for 
all  loaded  barges  and  the  number  of  coal  barges  is  currently 
assumed  to  be  a  constant  fraction  of  tow  size  but  these 
assumptions  are  easy  to  modify.  The  coal  barge  fractions 
vary  for  different  origin-destination  pairs.  Although  coal 
barge  fractions  are  constant  throughout  the  simulated 
period,  the  amount  delivered  from  each  tow  is  not  constant 
since  tow  sizes  are  randomly  generated,  as  mentioned 
previously. 

Cumulative  consumption  is  a  function  of  consumption 
rate  and  time.  The  mean  consumption  rate  is  constant  for 
each  utility  plant  during  each  simulation  period.  However, 
in  the  short  run  the  consumption  rate  fluctuates  randomly 
around  its  mean.  Currently,  the  consumption  rate  is  assumed 
to  be  uniformly  distributed.  The  actual  distribution  is  yet 
to  be  determined.  In  addition,  this  model  assumes  that  the 
consumption  rate  stays  constant  during  the  unit  time 
interval  which  is  specified  exogenously.  The  consumption 
rate  is  updated  every  time  unit  by  the  periodical  routine. 
Therefore,  cumulative  consumption  is  a  step-wise  linear 
function  over  time  whose  slopes  represent  consumption  rates. 

Inventory  levels  are  represented  in  this  model  by  the 
difference  between  cumulative  deliveries  and  cumulative 
consumption.  Whenever  inventory  levels  drop  to  negative 
values,  the  node  routine  computes  stock-out  amounts  and 
durations  for  the  analysis  of  total  costs. 

Periodical  Routine 


This  model  assumes  that  the  consumption  rate  stays 
constant  during  the  unit  time  interval  which  is  specified 
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exogenously.  The  periodical  routine  updates  the  consumption 
rate  every  time  unit,  and  also  updates  the  cumulative 
consumption,  the  inventory  status,  and  the  stock-out  amount 
and  duration  for  each  utility  plant  in  every  time  unit. 

There  are  three  functional  differences  between  the  node 
routine  and  the  periodical  routine. 

1.  The  periodical  routine  does  not  update  the 
cumulative  deliveries. 

2.  The  periodical  routine  updates  the  cumulative 
consumption  and  the  status  of  inventory  level  for 
each  utility  plant.  However,  the  node  routine 
updates  those  for  only  one  utility  plant  at  a 
time . 

3.  The  periodical  routine  can  update  the  consumption 
rate  based  on  a  seasonal  factor.  The  node  routine 
can  only  update  that  when  a  coal  tow  arrives  its 
destination. 

With  the  periodical  routine,  consumption  and  inventories  are 
estimated  more  precisely  than  without  it. 

Stall  Routine 


Stalls  are  random  failures  during  which  chambers  are 
not  available  to  serve  tows.  Thus  stalls  tend  to  increase 
delays  and  service  time  variability.  Stall  characteristics 
differ  among  chambers  and  are  defined  in  terms  of  durations 
and  frequencies  which  depend  on  weather  conditions  and 
physical  conditions  at  each  chamber.  Stalls  are  relatively 
rare  compared  to  other  events.  Their  occurrence  is  very 
difficult  to  predict.  Currently,  the  model  assumes  that  the 
time  intervals  between  two  successive  stalls  are 
exponentially  distributed. 

The  stall  routine  updates  the  available  chamber  time 
and  the  next  stall  occurrence  time  and  duration  at  this 
chamber.  The  available  chamber  time  is  defined  as  the  sum 
of  the  stall  occurrence  time  and  its  duration. 

Scheduler  Routine 


The  scheduler  routine  controls  the  simulation  clock  and 
invokes  the  five  operation  routines  to  process  the  necessary 
function  in  the  simulation  model.  The  procedures  in  the 
scheduler  routine  are  as  follows: 

1.  updating  the  next  event  occurrence  time  for  each 
event  type ; 
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2.  determining  the  next  event  occurrence  time  and 
event  type  among  the  five  types  of  events; 

3.  moving  the  simulation  clock  further  to  the 
occurrence  time  of  next  event;  and, 

4.  invoking  the  appropriate  operation  routine. 

Input  Requirements 

Generally,  four  types  of  inputs  are  required  to  operate 
the  simulation  model: 

1.  those  related  to  link  and  lock  characteristics; 

2.  those  related  to  traffic  demand  between  origins 

and  destinations; 

3.  those  related  to  probability  distributions;  and, 

4.  those  related  to  inventories  and  consumption. 

Link  and  Lock  Characteristics 

The  following  kinds  of  information  are  needed  for  each 

link; 

1 .  end  nodes ; 

2.  link  length; 

3.  distances  between  the  end  nodes  and  the  lock; 

4 .  number  of  chambers ; 

5.  average  frequencies  and  durations  of  stalls; 

6.  maximum  number  of  cuts  at  each  chamber; 

7.  average  service  time  for  each  number  of  cuts  at 
each  chamber; 

8.  maximum  number  of  barges  for  each  cut  size  at  each 
chamber; 

9.  bias  time  for  each  auxiliary  chamber;  and, 

10.  random  number  seeds  for  lock  service  times.  * 
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Traffic  Demand 


Traffic  demand  in  tows  per  day  is  expressed  in  the  form 
of  origin-destination  (0-D)  matrices  by  time  periods.  The 
lengths  of  time  periods  may  be  different  and  need  to  be 
specified.  The  required  input  data  for  traffic  demand  are 
as  follows: 

1.  origin  and  destination  nodes  for  each  0-D  pair; 

2.  average  trip  rate  in  tows  per  day  for  each  0-D 
pair,  each  direction  and  each  time  period; 

3 .  average  number  of  barges  per  tow  for  each  0-D 
pair  ; 

4.  fraction  of  coal  barges  in  a  tow  for  each  coal  0-D 
pair; 

5.  maximum  and  minimum  limits  of  tow  size  for  each 
coal  0-D  pair; 

6.  barge  payload  in  short -tons; 

7.  speed  distribution  (mean  and  standard  deviation); 

8.  ratio  of  backhaul  speed  to  linehaul  speed 
{full/empty  or  upstream/downstream) ;  and, 

9.  durations  of  time  periods. 

Probability  Distributions 

Probability  distributions. are  specified  in  this  model 
for  the  following: 

1.  lock  service  times 

2.  trip  generation 

3.  tow  composition  (barges  per  tow);  and, 

4.  coal  consumption  at  power  plants. 

The  probability  distribution  tables  represent 
cumulative  distribution  curves,  where  the  abscissas 
represent  cumulative  frequency,  and  the  ordinates  represent 
tfie  ratio  of  the  tabulated  variable  to  its  mean.  To  reduce 
the  input  com.plexity,  a  specified  number  of  equal  intervals 
is  currently  used  for  any  cumulative  frequency  distribution. 
Therefore,  only  the  values  on  the  ordinates  associated  with 
frequency  intervals  are  required. 
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Inventories  amd  Con«ii"ipHnn 


Initial  inventory  levels  in  short -tons  for  different 
nodes  (utility  plants)  must  be  specified.  In  addition, 
consumption  rates  in  short -tons  per  day  are  expressed  in  the 
form  of  node  matrices  by  time  period.  The  information  on 
cumulative  deliveries,  cumulative  consumption,  and  inventory 
levels,  is  provided  for  intervals  whose  duration  in  days 
must  be  specified. 

Model  Output 

This  model  provides  the  following  results: 

1.  mean  and  standard  deviation  of  waiting  time  at 
each  lock; 

2.  mean  and  standard  deviation  of  the  interarrival 
time  distribution  for  each  lock; 

3.  mean  and  standard  deviation  of  interarrival  time 
distribution  for  each  lock  and  direction; 

4 .  mean  and  standard  deviation  of  interdeparture  time 
distribution  for  each  lock; 

5.  mean  and  standard  deviation  of  interdeparture  time 
distribution  for  each  lock  and  direction; 

6.  total  tow  travel  time  (not  including  the  queuing 
time,  lockage  time,  and  dwell  time)  in  days; 

7.  mean  and  standard  deviation  of  travel  time  in 
hours  per  tow  between  any  lock  and  any  adjacent 
node ; 

8.  total  tow  travel  distances  in  1,000  miles; 

9.  average  travel  speed  in  miles  per  day; 

10.  average  lock  service  time  for  each  chamber  and 
cut  ; 

11.  total  number  of  tow  trips  for  each  0-D  pair; 

12.  monthly  cumulative  deliveries,  cumulative 
consumption,  and  inventory  levels  tables  in  1,000 
short-tons  for  different  power  plants; 

13.  cumulative  deliveries,  cumulative  consumption,  and 
inventory  levels  tables  for  specified  intervals  in 
1,000  short-tons  for  different  power  plants; 
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14.  graphs  of  cumulative  deliveries  and  cumulative 
consumption  by  specified  time  intervals  for 
different  power  plants;  and, 

15.  graphs  of  inventory  level  by  specified  time 
interval  for  different  power  plants. 


Validation 

To  check  the  logic  of  this  simulation  model,  its 
results  are  first  compared  to  theoretical  (but  very  well 
established)  results  from  queuing  theory.  This  also  checks 
the  ability  of  the  model  to  represent  general  stochastic 
distributions.  The  results  of  the  model  are  then  compared 
with  observed  data  to  demonstrate  how  closely  the  model 
represents  real  systems  and  verify  its  ability  to  simulate 
the  special  features  of  waterways. 

The  predicted  waiting  times  by  the  simulation  model  at 
a  single  lock  are  compared  with  those  obtained  from  queuing 
theory  when  arrivals  are  Poisson  distributed  and  service 
times  are  generally  distributed.  The  validation  is 
conducted  for  a  variety  of  volume/capacity  (V/C)  ratios 
ranging  from  0.0471  to  0.8934.  To  reduce  the  variance  of 
the  output  each  result  is  obtained  by  averaging  the  output 
from  30  independent  simulation  runs.  To  insure  results  are 
compared  for  a  steady  state,  each  simulation  run  discards 
the  first  10,000  tow  waiting  times  and  collects  the  next 
12,000  values  for  computing  the  average  waiting  time.  The 
results  are  shown  in  Table  1.  They  confirm  that  the 
simulated  and  theoretical  average  waiting  times  are 
extremely  close.  Such  results  verify  that  the  overall 
mechanism  of  the  simulation  model  is  correct.  They  also 
show  that  generally  distributed  service  times  are  generated 
satisfactorily  in  the  simulation  model.  That  is  reassuring 
since  the  same  logic  is  also  used  to  generate  generally 
distributed  interarrival  times  for  G/G/1  queues  and, 
ultimately  to  develop  metamodels  for  series  of  G/G/1  queues. 

The  simulation  results  are  then  compared  with  the 
observed  data  (from  January  1987)  at  Locks  22,  24,  25,  26, 
and  27  on  the  Mississippi  River.  These  particular  five 
locks  were  selected  mainly  because  they  were  considered 
especially  critical  to  the  entire  network  by  the  Corps  of 
Engineers.  Hence  extensive  data  analysis  and  performance 
evaluations  had  already  been  conducted  for  these  five  locks. 

At  that  time.  Locks  22,  24,  and  25  had  single  600  ft 
long  chambers.  Locks  26  and  27  had  two  chambers  (600  ft  and 
360  ft  long  at  Lock  26,  1200  ft  and  600  ft  at  Lock  27) .  It 
is  noted  that  the  simulation  results  for  the  five  lock 
system  were  obtained  simultaneously.  The  validation  results 
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TABLE  1.  COMPARISON  OF  THEORETICAL  AND  SIMULATED  RESULTS 


FOR  A  SINGLE 

LOCK 

QUEUE 

(M/G/1) 

v/c 

T 

■‘■A 

Avg  Var 

Avg 

T 

Var 

W  . 

Sim 

W^*'*  Devia. 

(hr)  (hr^) 

(hr) 

(hr^) 

(hr) 

(hr) 

(%) 

0.893 

0.888  0.789 

0.793 

0.319 

4 . 9516 

5 . 0059 

-1 . 09 

0 . 755 

0.888  0.789 

0.670 

0.227 

1 . 5575 

1.5522 

0 . 34 

0 . 566 

0.888  0.789 

0.503 

0.128 

0.4926 

0.4935 

-0.19 

0.330 

0.888  0.789 

0.293 

0.044 

0 . 1082 

0 . 1087 

-0.46 

0 . 047 

0.888  0.789 

0.042 

0.001 

0 . 00155 

0 . 00156 

-0.64 

*1  Ta 

*2  To 

*3  W  • 

"sim 

*4  W,. 

*5  Devia. 


interarrival  times 
service  times 

average  waiting  times  from  simulation 
average  waiting  times  from  queuing  theory 
deviation  which  is  defined  as  (Wg^^-W^) /W^*100% 


TABLE  2.  COMPARISON  OF  SIMULATED  AND  OBSERVED  AVERAGE 
WAITING  TIMES 


Lock 

w  . 

(min) 

W  K 
”obs 

(min) 

Difference 

(min) 

95%  Confidence 
Interval 

22 

4 . 09 

3.73 

0.36 

3.49 

24 

6 . 12 

6.36 

0.24 

6.72 

25 

4.49 

10.94 

6.45 

-  *3 

26 

119.40 

130.99 

11.59 

60.73 

27 

36.49 

34.43 

2.06 

23 . 92 

*1  Wgin,  :  simulated  average  waiting  times 
*2  :  observed  average  waiting  times 

*3  The  comparison  is  not  appropriate. 
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TABLE  3.  COMPARISON  OF  CHAMBER  VOLUMES 


Lock 

Chamber 

Volsi.'" 

(tows/day) 

Volobs*" 

(tows/day) 

Difference 

(tows/day) 

95%C.  I  . 

22 

1 

44.61 

45 

0 .39 

1 . 88 

24 

1 

56 . 00 

56 

0 . 00 

2.24 

25 

1 

51.40 

52 

0 . 60 

2 . 13 

26 

1 

275.20 

265 

10.20 

2 . 83 

26 

4 

155.96 

167 

11.04 

3 . 81 

27 

1 

390.95 

389 

1.95 

10.25 

27 

4 

306.73 

306 

0.7 

10.62 

*1  Volgip,  :  simulated  volumes 
*2  VoIqJjq  :  observed  volumes 

*3  95%C.I.:  95%  confidence  interval  based  on  t  test 


TABLE  4.  COMPARISON  OF  CUT  VOLUMES 


Lock 

Cham*^ 

Cuts 

Vol  • 

(tows/day) 

(tows/day) 

Difference 

(tows/day) 

95% 

C.I. 

22 

1 

1 

30.35 

31 

0.65 

1.47 

22 

1 

a2 

14.26 

14 

0.26 

- 

24 

1 

1 

39.76 

40 

0.24 

1 . 76 

24 

1 

&2 

16.24 

16 

0.24 

- 

25 

1 

1 

35.88 

37 

1.12 

1 . 65 

25 

1 

s2 

15 . 52 

15 

0 . 52 

0 . 97 

26 

1 

1 

75.46 

74 

1.46 

2 . 04 

26 

1 

^2 

199.74 

191 

8 . 74 

3 .25 

26 

4 

1 

137.25 

147 

9 . 75 

2 . 87 

26 

4 

fe2 

18 . 71 

20 

1.29 

- 

27 

1 

1 

390 . 95 

389 

1 . 95 

10.25 

27 

4 

1 

269 . 05 

265 

4 . 05 

5 . 59 

27 

4 

z2 

37.67 

41 

3.33 

5 . 76 

*1  Cham  :  chamber 

*2  "Cuts"  :  are  subsets  of  barges  into  which  tows  are 

subdivided  for  passage  through  lock  chambers 
*3  Volgj^n,  :  simulated  volumes 

*4  Volgi^g  :  observed  volumes 

*5  95%  C.I.  ;  95%  confidence  interval  based  on  t  test. 
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are  summarized  in  Tables  2,  3,  and  4.  Each  result  is 
averaged  from  80  independent  simulation  runs.  The  initial 
condition  for  simulation  is  assumed  to  be  an  empty 
system, which  is  consistent  with  the  observed  condition  for 
this  system  in  winter. 

Table  2  shows  that  the  difference  between  the  simulated 
and  observed  average  waiting  times  for  each  lock  are  within 
the  95%  confidence  interval  based  on  the  t  test,  except  at 
Lock  25.  The  observed  data  also  show  that  tows  sometimes 
were  kept  waiting  at  Lock  25  even  when  the  chamber  was  idle. 
Such  operation  is  somewhat  unusual.  Therefore,  no  direct 
comparison  of  average  waiting  times  at  Lock  25  is 
appropriate.  Tables  3  and  4  also  show  that  the  simulation 
model  represents  the  real  system  quite  well. 

Each  simulation  run  takes  from  a  few  seconds  to  a  few 
minutes  on  a  personal  computer,  depending  on  traffic 
volumes,  duration  of  simulation  periods,  network  size,  and 
other  factors.  Despite  that,  simulation  time  becomes 
expensive  for  evaluating  large  combinatorial  investment 
scheduling  problems.  For  example,  when  there  are  n=20 
possible  investment  projects,  it  is  necessary  to  simulate 
2^°  combinations  to  make  the  best  decision.  30*2^°  separate 
simulation  runs  are  then  required  if  each  performance 
measure  is  based  on  the  average  over  30  independent 
replications.  Furthermore,  the  project  combinations  may 
have  to  be  evaluated  over  several  time  periods.  Therefore, 
as  n  increases,  direct  evaluation  by  simulation  becomes  very 
expensive . 

A  metamodeling  approach  is  proposed  to  overcome  the 
computational  requirements  of  simulation.  A  simulation 
model  can  then  be  treated  as  a  function  with  unknown 
explicit  form  that  turns  input  parameters  into  output 
performance  measures.  The  metamodeling  approach  provides  a 
method  to  develop  simple  formulas  to  approximate  this 
function. 


35 


36 


CHAPTER  4 
NUMERICAL  METHOD 


Methodology 

In  this  study,  a  numerical  method  has  been  developed 
for  estimating  delays  through  a  series  of  queues.  This 
method  was  originally  developed  for  systems  with  bi¬ 
directional  servers.  With  a  few  simplifications,  this 
method  can  be  adapted  for  the  more  generally  encountered 
systems  with  one-directional  servers. 

The  numerical  method  decomposes  large  systems  into 
single  lock  queuing  stations  and  then  analyzes  the 
interarrival  and  interdeparture-time  distributions  at  theses 
single  locks  one  by  one.  Although  the  numerical  method  is  a 
decomposition  procedure,  each  lock  is  not  analyzed 
independently.  The  interarrival-time  distribution  at  one 
lock  is  affected  by  the  interdeparture -time  distributions  at 
ad j  acent  locks . 

The  method  consists  of  three  major  modules:  arrival 
processes,  departure  processes,  and  delay  functions. 

Arrival  processes  at  a  particular  lock  depend  on  the 
interdeparture  time  distributions  from  the  upstream  and 
downstream  locks.  The  departure  processes  depend  on  the 
interactions  among  the  interarrival -time  distributions  and 
service- time  distributions  at  one  lock.  The  delay  functions 
define  the  relations  among  waiting  times,  interarrival -time 
distributions,  interdeparture -time  distributions  and 
service-time  distributions.  The  basic  concept  of  this 
method  is  to  identify  two  parameters,  namely  the  mean  and 
coefficient  of  variation,  of  the  interarrival  and 
interdeparture-time  distributions  for  each  lock,  and  then 
estimate  the  implied  waiting  times.  Currently,  the 
following  assumptions  are  used  in  the  numerical  method. 

1.  Interarrival  times  and  service  times  are  generally 
(i.e.,  arbitrarily)  distributed. 

2.  Each  lock  has  one  chamber. 

3  .  Inflows  and  outflows  occur  only  at  the  two  end- 
nodes  of  a  series  of  locks. 

4 .  The  average  upstream  volumes  are  equal  to  the 
downstream  volumes. 

5.  The  volume-to-capacity  ratio  (p)  is  less  than  1.0 
at  every  lock. 


37 


It  should  be  noted  that  Assumptions  2,  3,  4,  and  5  are 
only  applicable  to  the  numerical  method.  The  simulation  model 
is  not  limited  by  those  assumptions.  The  numerical  method  can 
provide  a  quick  and  inexpensive  approach  for  the  analysis  of 
lock  delays.  However,  Assumptions  2,  3,  4,  and  5  limit  the 
applicability  of  the  currently  developed  numerical  method  and 
necessitate  the  substitution  of  the  simulation  model  when 
significant  deviations  from  those  assumptions  must  be 
considered.  With  some  extensions  to  the  numerical  method 
expected  in  the  near  future,  Assumptions  2  and  3  may  be 
eliminated.  Assumption  4  cculd  be  relaxed  even  though  it  is 
usually  realistic  for  waterways. 

Each  module  of  the  numerical  method  consists  of  one  or 
more  metamodels  that  are  functional  relations  whose  parameters 
are  statistically  estimated  from  simulation  results.  Thus, 
simulation  is  not  just  used  to  validate  these  metamodels;  it 
is  the  basis  for  their  development. 

Two  simulation  experiments  were  conducted  in  this  study. 
Experiment  1  served  two  goals  here,  namely  the  development  of 
metamodels  and  the  validation  of  these  metamodels.  A  split- 
sample  analysis  was  conducted  to  achieve  both  of  the  above 
goals  in  a  single  experiment.  The  split-sample  analysis 
proceeded  as  follows.  First,  the  collected  data  points  from 
Experiment  1  were  grouped  in  strata  based  on  one  or  more 
important  categorical  variables  (This  method  is  called  a 
stratum-specific  split-sample  scheme  [20] ) .  Second,  all  data 
points  within  a  stratum  were  randomly  assigned  to  one  of  two 
groups,  namely  the  training  group  (used  to  develop  metamodels) 
and  the  holdout  group  (used  to  validate  metamodels  and  to  test 
their  reliability) .  Such  assignment  was  done  separately  for 
each  stratum.  The  goal  of  stratum-specific  random  assignment 
is  to  insure  that  the  two  groups  of  data  (training  and 
holdout)  are  equally  representative  of  the  parent  population 
[20]  .  In  this  study,  a  computer  randomly  assigned  2/3  of  data 
to  the  training  group  and  1/3  of  data  to  the  holdout  group. 
It  is  noted  that  the  above  two  steps  were  conducted  before  any 
data  were  analyzed. 

Experiment  2  aimed  to  validate  the  numerical  method 
against  the  simulation  model.  The  ranges  of  variables  in 
these  experiments  were  carefully  chosen  to  reflect  the  ranges 
encountered  in  waterway  systems. 

The  procedures  used  in  developing  each  metamodel  are 
summarized  below. 

1.  Queuing  theory  was  applied  to  identify  the  input 
(independent)  variables  which  will  affect  the  output 
measures  (dependent  variables) ,  and  to  propose 
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variables  and  the  output  measures  with  appropriate 
structure  form. 

2.  The  relevant  output  measures  (dependent  variables) 
were  plotted  versus  the  input  (independent) 
variables.  This  step  helped  confirm  the  relations 
between  the  output  measures  and  the  input 
parameters . 

3.  Pearson  correlation  coefficients  were  computed. 

This  step  confirmed  correlations  between  selected 
dependent  variables  and  the  independent  variables 
and  helped  avoid  multicollinearity  problems  among 
the  independent  variables. 

4 .  The  parameters  were  estimated  for  the  proposed 
functional  relations  ("metamodels").  The  selection 
of  the  preferred  metamodel  was  based  mainly  on  the 
sample  squared  multiple  correlation  (R^) ,  the 
standard  error  of  residuals,  and  the  residual 
analysis.  The  defines  the  explanatory  power  of 
the  alternative  metamodels.  In  general,  the 
metamodel  with  R^  closest  to  1  is  preferred  since 
it  best  accounts  for  the  variation  of  dependent 
variables.  It  is  also  important  to  examine  if  the 
independent  variable  is  significant  in  explaining 
the  variation  of  dependent  variable. 

5.  Residual  analysis  was  performed  to  detect  outliers 
and  to  check  whether  any  metamodel  violated  certain 
regression  assumptions,  such  as  normality  and 
homoscedasticity .  The  residual  analysis  should 
include  the  property-analysis  and  the  graphical 
analysis  of  residuals.  The  basic  residual 
properties  to  be  examined  include  the  mean, 
variance,  skewness,  and  kurtosis.  The  residual 
mean  should  equal  0.  The  variance  is  the  residual 
mean  square .  Skewness  indicates  the  degree  of 
asymmetry  of  a  distribution;  it  should  be  close  to 
0  to  avoid  violating  the  normality  assumption. 
Kurtosis  indicates  the  heaviness  of  the  tails 
relative  to  the  middle  of  the  distribution;  it 
should  be  close  to  3  to  avoid  violating  the 
normality  assumption.  (The  commercial  statistical 
software  usually  deducts  3  from  the  value  of 
kurtosis.  Therefore,  the  kurtosis  should  be  close 
to  0  if  it  is  obtained  from  commercial  statistical 
software.)  However,  skewness  and  kurtosis 
statistics  are  highly  variable  in  small  samples  and 
hence  are  often  difficult  to  interpret  [20]  .  The 
graphical  analysis  is  the  most  direct  and  revealing 
way  to  examine  a  set  of  residuals.  The  residuals 
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can  be  displayed  in  one  or  two  dimensions.  The 
useful  one-dimensional  plot  of  residuals  includes 
histograms  (or  stem-and-leaf ) ,  schematic  plots,  and 
normal  probability  plots  [20] .  The  two-dimensional 
plot  examines  the  relationships  of  the  residuals  to 
either  dependent  or  independent  variables  and  is 
sometimes  useful  for  identifying  violation  of 
regression  assumptions  [20] . 

Structure  c£  Numerical  Method 

Basically,  the  numerical  method  is  a  decomposition 
model .  It  decomposes  systems  of  queues  into  separate 
queuing  stations.  The  analysis  of  each  queuing  station  is 
decomposed  into  three  steps,  namely  arrival  processes, 
departure  processes  and  delays. 

ihe  arrival  processes  module  serves  two  functions. 
First,  for  each  direction  at  each  lock  it  identifies  the 
interarrival -time  distribution  based  on  the  interdeparture¬ 
time  distribution  from  the  previous  lock  and  the  intervening 
speed  distribution.  Second,  the  overall  interarrival -time 
distribution  at  each  lock  is  estimated  by  combining  the 
interarrival-time  distributions  from  the  two  adjacent  locks. 
The  arrival  processes  module  is  very  important  in 
identifying  and  taking  into  account  the  interdependence 
among  locks,  even  when  the  system  is  decomposed  into 
individual  lock  queuing  stations. 

The  departure  processes  module  also  serves  two 
functions.  First,  the  overall  interdeparture-time 
distribution  at  each  lock  is  estimated  based  on  the  overall 
interarrival- time  distribution  and  the  lock's  service-time 
distribution.  Second,  the  overall  interdeparture- tine 
distribution  at  each  lock  is  split  into  directional 
interdeparture-time  distributions.  Thus,  at  each  lock,  the 
departure  processes  module  estimates  two  output 
distributions  (interdeparture  times  for  two  directions) 
based  on  three  input  distributions  (interarrival  times  for 
two  directions  and  bi-directional  service  times) . 

The  delay  module  estimates  average  waiting  time  at  each 
lock.  The  inputs  for  the  delay  function  are  the  variances 
of  the  interarrival  times,  interdeparture  times,  and  service 
times.  While  the  service-time  distributions  at  each  lock 
are  exogenously  specified  inputs,  the  parameters  of  the 
interarrival-time  and  interdeparture-time  distributions  are 
obtained  from  the  arrival  processes  module  and  the  departure 
processes  module,  respectively. 

A  scanning  procedure  is  employed  in  the  numerical 
method  to  analyze  successive  locks.  Interrelations  among 
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the  locks  in  a  system  are  accounted  for  since  the  departure 
distribution  from  one  lock  influx  aces  the  arrival 
distribution  at  the  next  lock.  In  the  initial  scan  along 
one  direction,  no  endogenous  estimate  is  yet  available  of 
arrival  distributions  from  the  opposing  direction  (unless 
there  is  only  one-way  traffic) .  If  there  is  opposing 
traffic,  some  initializing  estimates  must  be  provided  for 
the  variance  of  interdeparture  times  from  the  opposing 
direction.  The  scans  are  then  repeated  in  an  iterative 
algorithm  which  alternates  scanning  direction.  After  the 
first  iteration,  the  initializing  estimates  are  replaced 
with  endogenous  estimates  from  previous  iterations.  The 
iterative  scans  continue  until  convergence  is  achieved  in 
values  of  variables  such  as  wait  times  or  interdeparture - 
time  variances. 

Sampling  Plan 

The  data  points  in  Experiment  1  were  collected  from  40 
different  simulated  lock  systems.  Each  system  consisted  of 
three  locks.  The  trip  rate  differed  for  each  lock  system 
and  remained  constant  throughout  the  data  collection  period. 
Therefore,  the  data  points  in  Experiment  1  represent  the 
characteristics  of  40  different  trip  rates.  These  40 
different  trip  rates  were  generated  uniformly  in  the  range 
between  0.0417  tows/hr  (1  tow/day)  and  1.6667  tows/hr  (40 
tows/day) .  The  range  was  selected  to  cover  all  possible 
levels  empirically  found  on  waterways. 

The  three  locks  within  any  system  differed  from  each 
other,  in  terms  of  their  Volume  to  Capacity  (V/C)  ratio  and 
the  squared  coefficient  of  variation  of  service  times. 
Therefore,  the  data  points  in  Experiment  1  represent  the 
characteristics  of  120  V/C  ratios  and  120  squared 
coefficients  of  variation  of  lock  service  times.  The  V/C 
ratios  were  generated  uniformly  between  0.04  to  0.97, 
covering  all  observed  waterway  operations.  The  squared 
coefficients  of  variation  of  lock  service  times  were 
generated  uniformly  between  0.34  to  0.81,  covering  almost 
all  waterway  operations. 

The  distances  between  locks  and  the  tow  speed 
distributions  were  also  different  in  each  lock  system.  The 
distances  were  generated  uniformly  between  0  to  156  miles. 
The  average  tow  speeds  were  uniformly  distributed  between  72 
to  295  miles/day.  The  standard  deviations  of  tow  speeds 
were  uniformly  distributed  between  2.22  to  136.3  miles/day. 
The  ranges  of  distances  and  speed  distributions  were 
selected  to  cover  most  waterway  situations  of  practical 
interest . 
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At  each  lock,  the  collected  characteristics  for  its 
corresponding  data  point  include  the  directional 
interarrival -time  distributions,  the  directional 
interdeparture -time  distributions,  the  overall  interarrival - 
time  distribution,  the  overall  interdeparture-time 
distribution,  and  the  average  waiting  time.  The 
corresponding  trip  rate,  V/C  ratio,  lock  service  time 
distribution,  distances,  and  tow  speed  distribution  were 
also  recorded.  Therefore,  the  database  ii.  Experiment  1 
contains  120  different  average  waiting  times,  120  different 
overall  interarrival -time  distributions,  120  different 
overall  interdeparture-time  distributions,  240  different 
directional  interarrival-time  distributions,  and  240 
different  directional  interdeparture- time  distributions. 

To  obtain  the  above  data,  each  lock  system  was 
simulated  with  30  independent  replications.  To  ensure  that 
the  results  were  for  a  steady  state,  each  simulation  run 
discarded  the  first  10,000  observations  and  collected  the 
next  12,000  values  for  computing  the  data  characteristics. 
10,000  observations  were  discarded  based  on  moving  average 
analysis  [33] .  In  addition,  the  random  number  seeds  were 
all  200,000  apart,  which  ensured  the  random  number  streams 
would  be  independent  since  each  independent  replication 
discarded  10,000  observations  and  collected  the  next  12,000 
values.  Therefore,  the  average  waiting  time  at  each  lock 
was  obtained  by  the  following  steps. 

1.  Thirty  independent  replications  of  the  simulation 
were  made  for  each  lock  system,  each  generating  an 
output  sequence  of  22,000  tows.  Hence,  a  set  of 
observations  were  obtained: 

Wj^=  n-th  member  of  the  output  (tow  waiting 

time)  sequence  from  the  m-th  independent 
replication  of  the  simulation,  m=l, . . . 
30,  n=l, . . .  22,000 . 

2.  The  sample  mean  of  the  average  waiting  time  was 
obtained  from  each  replication  [33]  . 


22,000 


K  = 


12, 000 


n=10, 001 


(9) 


3 .  The  unbiased  estimate  of  the  average  waiting  time  W 
is  estimated  as  follows  [33] : 
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(10) 


1 

30 


30 


/n=l 


The  variances  of  overall  interarrival -time  distributions 
were  obtained  by  the  following  five  steps. 

1.  A  set  of  observations  were  obtained  from  30 
independent  simulation  replications: 


tAmn=  n-th  member  of  the  overall 

interarrival -time  sequence  from  the 
m-th  independent  replication  of  the 
simulation,  m=l, . . .30, 
n=l, . . .22, 000 . 


2.  The  sample  mean  of  the  overall  interarrival  times 
was  obtained  on  each  replication  [33]  : 


t 


Am 


1 

12,000 


22,000 


E 

71=10,001 


'Amn 


(11) 


3  . 


The  estimate  of  the  average  interarrival  time  was 
obtained  as  follows  [33] : 


30 


1  -s  o  ^Am 
in=l 


(12) 


4  . 


The  variance  estimate  was  obtained  for  each 
replication  [33] ; 


22,000 


^An\ 


12/000  ^=10,001 


^  ^Amn 


(13) 


5.  The  variance  of  overall  interarrival  times  was 
obtained  as  follows  [33] : 


(14) 


The  variances  of  overall  interdeparture- time  distributions, 
directional  interarrival-time  distributions,  and  directional 
interdeparture- time  distributions  also  were  obcained 
following  the  above  procedures. 
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Arrival  Procasses 


The  following  two  steps  were  used  for  estimating  the 
mean  and  standard  deviations  of  interarrival  times. 


Step  1 

Estimate  the  means  and  standard  deviations  of 
directional  interarrival  times  at  a  particular  lock  from  the 
directional  interdeparture -time  distributions  of  the 
adjacent  upstream  and  downstream  locks. 


If  flows  are  conserved  between  locks  and  if  the  V/C 
ratio  is  less  than  1,  the  average  directional  arrival  rates 
at  one  lock  should  be  equal  to  the  average  directional 
departure  rates  from  adjacent  upstream  and  downstream  locks. 
Therefore,  the  average  directional  interarrival  times  at 
that  lock  should  also  be  equal  to  the  average  directional 
interdeparture  times  from  adjacent  upstream  and  downstream 
locks.  Such  relations  are  represented  by: 


ik=i-l,  if  j=l 
*ic=i+l,  if  j=2 


where 

7  :  the  average  interarrival  time  for  Direction  j  and 

Lock  i 

7  :  the  average  interdeparture  time  for  Direction  j  and 

Lock  k 

j  :  direction  index  (1  =  downstream,  2  =  upstream) 

If  each  tow  moves  at  the  same  speed,  the  directional 
interarrival  time  distributions  at  one  lock  will  be 
identical  to  the  directional  interdeparture  time 
distributions  at  the  preceding  lock.  However,  speed 
variations  change  headway  distributions  along  the  distance 
between  locks. 

A  metamodel  is  developed  to  estimate  the  standard 
deviation  of  directional  interarrival  times  at  one  lock. 
Before  developing  the  metamodel,  the  database  obtained  in 
Experiment  1  was  split  into  two  groups  of  data  (training  and 
holdout) .  The  stratum-specific  split-sample  scheme  was 
employed  for  data  assignment.  For  developing  the  necamodel, 
4  data  points  were  available  from  each  of  the  40  iock 
systems  (downstream:  the  relation  between  Locks  1  and  2,  and 
the  relation  between  Locks  2  and  3;  upstream:  the  relation 
between  Locks  3  and  2,  and  the  relation  between  Locks  2  and 
1) .  Therefore,  a  total  of  160  data  points  was  available 
from  Experiment  1  for  metamodel  development.  The  split 
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ratio  between  the  training  and  holdout  groups  is 
approximately  2:1.  Therefore,  107  observations  were 
included  in  the  training  group  and  53  observations  were 
included  in  the  holdout  group. 


The  following  metamodel  was  estimated  with  data  from 
the  training  group  is  as  follows: 

^iiPviks  ik=i-l,if  j  =  l 

'Jc=i+l,if  j=2  (16) 

(0.002) 


aji 


--a^j^  +  0 . 02511n  (1  +  - 


=  0.999954 


n  =  107 


Sg  =  0.0586 


=  5.168561 


where 


standard  deviation  of  interarrival  times  for 
Direction  j  and  Lock  i 


:  standard  deviation  of  interdeparture  times  for 

Direction  j  and  Lock  k 


£)..  :  distance  between  Locks  i  and  k 

ik 

M’vi/c  •  average  tow  speed  between  Locks  i  and  k 


^vik 

j 


standard  deviation  of  tow  speeds  between  Locks  i 
and  k 

direction  index  (1  =  downstream,  2  =  upstream) 


S 


e 


Standard  error  of  dependent  variable 


My 


mean  of  dependent  variable 


The  numbers  shown  in  the  parentheses  are  the  standard  errors 
of  the  estimated  parameters  directly  above. 


Currently,  there  is  less  theoretical  basis  for  this 
metamodel  than  for  the  other  metamodels  developed  in  this 
study.  This  metamodel  was  developed  largely  by  empirical 
analysis.  The  dependent  variable  (standard  deviation  of 
directional  interarrival -time  distribution)  was  plotted 
versus  possible  influential  factors  that  include  the 
standard  deviation  of  directional  interdeparture -time 
distribution,  distance  between  two  locks,  average  tow  speed, 
and  standard  deviation  of  tow  speeds.  The  correlation 
coefficients  between  the  dependent  variable  and  influential 
factors  were  also  computed.  The  plot  and  the  correlation 
coefficient  show  that  the  standard  deviation  of  the 
directional  interarrival -time  distribution  is  highly 
correlated  to  the  standard  deviation  of  the  directional 
interdepa:' ture-time  distribution.  However,  the  plot  and  the 
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correlation  coefficient  also  show  that  the  standard 
deviation  of  the  directional  interarrival -time  distribution 
is  also  affected  by  the  factors  of  distance,  average  tow 
speed,  and  standard  deviation  of  tow  speeds.  Therefore, 
various  structural  forms  were  considered  for  this  metamodel, 
and  two  were  intensively  pursued. 

The  first  metamodel  for  includes  only  one 
independent  variable  (the  standard  deviation  of  the 
directional  interdeparture -time  distribution) .  The  second 
metamodel  for  also  considers  the  distance,  average  tow 
speed,  and  standard  deviation  of  tow  speeds  as  influential 
factors.  The  for  the  first  metamodel  is  0.999885.  The 
descriptive  statistics  and  test  results  for  the  first 
metamodel  are  included  in  Table  5. 

The  structural  form  of  the  second  metamodel  was 
preferred  since  it  satisfies  an  important  logical 
constraint.  The  standard  deviation  of  the  directional 
interarrival  times  should  be  equal  to  the  standard  deviation 
of  the  directional  interdeparture  times  from  the  preceding 
lock  when  the  distance  between  these  two  locks  is  0  or  when 
the  standard  deviation  of  tow  speeds  is  0.  The  chosen 
second  type  metamodel  (Eq.  16)  does  satisfy  those  logical 
constraints . 

It  is  noteworthy  that  in  Eq.  16,  the  coefficient  for 
the  standard  deviation  of  directional  interdeparture  times 
from  the  preceding  lock  is  equal  to  1  and  the  intercept  is 
0.  In  fact,  in  the  first  version  of  this  metamodel,  the 
intercept  was  0.00526  with  a  standard  error  of  0.021  which 
indicated  that  the  intercept  was  not  significant.  The 
coefficient  for  the  standard  deviation  of  directional 
interdeparture  times  was  0.999,  or  approximately  equal  to  1. 
This  suggested  that,  theoretically,  the  standard  deviation 
of  directional  interarrival  times  should  be  equal  to  the 
standard  deviation  of  directional  interdeparture  times  plus 
an  adjustment  factor  depending  on  the  speed  distribution  and 
distance.  Therefore,  this  metamodel  was  redeveloped  under 
two  restrictions:  (1)  the  intercept  should  be  0  and  (2)  the 
coefficient  for  the  standard  deviation  of  directional 
interdeparture  times  should  be  1.  The  descriptive 
statistics  and  test  results  for  the  second  type  metamodel 
are  listed  in  Table  6 . 

The  comparison  between  Tables  5  and  6  shows  that  the 
second  meta  del  does  predict  better  the  standard  deviation 
of  directional  interarrival- time  distribution.  The  of 
the  second  metamodel  is  slightly  higher.  However,  the 
statistics  of  £  and  £%  of  the  second  metamodel  are 
much  better.  The  mean  of  £  is  0.04  for  the  second  and  0.07 
for  the  first.  The  mean  of  £%  is  1.68  for  the  second  and 
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TABLE  5.  DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
TRAINING  GROUP,  FIRST  METAMODEL  FOR 


Variable 

Mean 

S.D.* 

Min 

Max 

""aji 

5 . 17 

8 . 61 

1.17 

49 . 57 

5.10 

8.67 

1.03 

49 . 57 

£ 

0 . 07 

0.07 

0.00041 

0 .29 

£%*** 

3.17 

3 . 02 

0.0016 

12.66 

n-"  =  107, 

0 . 999885 

*S.D.  :  standard  deviation 

**£  :  absolute  error  of  between  simulation  and 

metamodel 

***£%  :  percentage  error  of  between  simulation  and 

metamodel 

+n  :  number  of  observations 

++R]^^  :  the  coefficient  of  determination 


TABLE  6.  DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
TRAINING  GROUP,  SECOND  METAMODEL  FOR 

at 


Variable 

Mean 

S.D.* 

Min 

Max 

^aji 

5 . 17 

8.61 

1.17 

49 . 57 

‘^djk 

5 . 10 

8.67 

1.03 

49.57 

Dik 

77.83 

46.65 

0.00 

156.00 

^vik 

16:^. 00 

56.25 

72.00 

295 . 00 

47.54 

32 .30 

2.22 

136.30 

£ 

0 . 04 

0.04 

0.00046 

0.21 

£%*** 

1.68 

1.57 

0.0016 

7.60 

n-"  =  107, 

=  0.999954 

*S.D.  :  Standard  deviation 

**£  :  absolute  error  of  between  simulation  and 

metamodel 

***£%  :  percentage  error  of  between  simulation  and 

metamodel 

+n  :  number  of  observations 

+  +R]^^  :  the  coefficient  of  determination 
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3.17  for  the  first.  The  maximum  £  is  0.21  for  the  second 
and  0.29  for  the  first.  The  maximum  £%  is  7.6  for  the 
second  and  12.66  for  the  first.  Therefore,  the  second 
metamodel  was  chosen  for  predicting  the  standard  deviation 
of  the  directional  interarrival  times. 


The  holdout  group  data  was  used  to  test  the  reliability 
of  this  metamodel,  which  is  used  to  predict  the  standard 
deviation  of  directional  interarrival  times.  It  is  assumed 
in  this  test  that  the  values  obtained  from  the  simulation 
model  are  exact.  Therefore,  the  deviations  between  the 
values  obtained  from  simulation  and  metamodel  are  considered 
errors.  The  percentage  error  is  defined  as  the  error 
divided  by  the  value  obtained  from  simulation.  In  this 
test,  the  absolute  values  are  used  to  compute  errors  and 
percentage  errors.  The  use  of  this  metamodel  is  to  predict 
the  standard  deviation  of  directional  interarrival  times  for 
the  holdout  group  constitutes  cross-validation.  The 
descriptive  statistics  and  test  results  for  the  holdout 
group  are  listed  in  Table  7. 

The  shrinkage  on  cross-validation,  defined  as  - 
R2^,  indicates  the  reliability  of  metamodels.  In  general, 
shrinkage  values  less  than  0.10  are  indicative  of  a  reliable 
model  [20] .  Therefore  this  metamodel  is  reliable  since  the 
shrinkage  is  0.000001.  In  addition,  the  absolute  error 
(mean=0.04  and  maximum  value=0.17)  and  absolute  percentage 
error  (mean=1.74  and  maximum  value=7.76)  in  Table  7  also 
show  this  metamodel  predicts  the  standard  deviation  of 
directional  interarrival  times  well.  However,  the  maximum 
£%  is  7.6  for  the  training  group  and  7.76  for  the  holdout 
group,  which  indicates  that  the  directional  arrivals 
metamodel  sometimes  generates  considerable  errors  and  leaves 
room  for  improvement . 

Step  2 


The  mean  and  coefficient  of  variation  of  the  overall 
interarrival -time  distribution  for  this  lock  are  estimated 
baccd  on  the  coef  f  i  cici.t  ^.^f  variation  of  directional 
interarrival  times.  The  coefficient  of  variation  of 
directional  interarrival  times  could  be  obtained  by  dividing 
the  standard  deviation  of  directional  interarrival -time 
distribution  with  its  mean. 


t 


Ai 


(17) 


=  0.179+0.41(0,,/  +  0,,/) 
(0.027)  fn.oi4) 

R^  =  0.9188  n  =  79  s^  =  0.0059 


(18) 

pLy  =  0.988 
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the  average  interarrival  time  at  Lock  i 


squared  coefficient  of  variation  of  interarrival 
times  at  Lock  i 

squared  coefficient  of  variation  of  directional 
interarrival  times  for  Direction  j  and  Lock  i 

The  meaning  of  Eq.  17  would  be  clearer  if  viewed  in 
terms  of  the  average  trip  rates  rather  than  the  average 
interarrival  times.  Eq.  17  implies  that  the  overall  arrival 
rate  at  certain  lock  is  the  sum  of  the  average  directional 
arrival  rates  from  upstream  and  downstream. 

Before  developing  the  metamodel  of  Eq.  18,  the  stratum- 
specific  split-sample  scheme  was  also  employed  to  split  the 
database  of  Experiment  1  into  the  two  groups  (training  and 
holdout) .  For  developing  the  metamodel,  there  were  120 
observations  available  from  the  database  of  Experiment  1. 
Each  observation  was  obtained  from  one  individual  lock  and 
included  the  characteristics  of  overall  interarrival -time 
distribution  and  corresponding  directional  interarrival -time 
distributions  at  one  lock.  The  split  ratio  between  the 
training  and  holdout  groups  was  approximately  2:1. 

Therefore,  the  computer  randomly  assigned  79  observations  to 
the  training  group  and  41  observations  to  the  holdout  group. 

In  Eq.  18,  the  squared  coefficients  of  variations  of 
upstream  and  downstream  interarrival  times  carry  the  same 
weight  in  estimating  the  overall  variance  of  interarrival 
times,  since  directional  trip  rates  are  equal  according  to 
Assumption  4.  Eq.  18  should  be  reestimated  when  applied  to 
a  more  general  network  of  queues  with  imbalanced  flows. 

The  data  in  the  holdout  group  was  also  used  to  validate 
the  reliability  of  this  metamodel  and  the  cross-validation 
test  was  again  conducted.  The  descriptive  statistics  and 
test  results  for  the  training  group  and  the  holdout  group 
are  listed  in  Tables  8  and  9. 

The  shrinkage  (Rj^^-R2^)  in  this  metamodel  is  0.0282 
which  indicates  that  this  metamodel  is  fairly  reliable. 

This  becomes  more  convincing  after  examining  the  absolute 
error  and  absolute  percentage  error  in  the  cross-validation 
test.  The  absolute  error  ranges  from  0.00022  to  0.03  with  a 
mean  of  0.007.  The  absolute  percentage  error  ranges  from 
0.022  to  3.26  with  a  mean  of  0.70. 


where 


^Ai 


C  ^ 

'-Al 


C  ■  ^ 
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TABLE  7.  DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
HOLDOUT  GROUP  (a^^) 


Variable 

Mean 

S.D.* 

Min 

Max 

‘^aji 

5 . 04 

7.82 

1.23 

49.10 

*^81  k 

4.98 

7.84 

1.08 

49.10 

^ik 

78.34 

45.64 

0 . 00 

152.00 

/^vik 

171.56 

60 . 09 

79.20 

295.00 

^vik 

45.78 

34.54 

2.22 

136.30 

0.04 

0.04 

0 . 00004 

0 . 17 

£%*** 

1.74 

1.81 

0.0027 

7.76 

n-"  =  53, 

"D  2  +  +  _ 

_ji2 _ L. 

0.999953 

*S.D.  :  standard  deviation 

**£  :  absolute  error  of  between  simulation  and 

metamodel 

***£%  ;  percentage  error  of  between  simulation  and 

metamodel 

+n  :  number  of  observations 

++R2^  :  the  coefficient  of  determination 


TABLE  8.  DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
TRAINING  GROUP  (C;^^) 


I  Variable 

Mean 

S.D, 

Min 

Max 

c  •  ^ 

0 . 99 

0.02 

0 . 89 

1.02 

r  .  ^ 

0 . 98 

0 . 04 

0 . 74 

1.02 

0.99 

0.03 

0.84 

1 . 02 

0 . 004 

0 . 004 

0.00007 

0 . 02 

0.46 

0.40 

0.0066 

1.79 

n^  =  79, 

= 

0 . 9188 

*S.D.  :  standard  deviation 

**£  :  absolute  error  of  between  simulation  and 

metamodel 

***£%  :  percentage  error  of  CAi"  between  simulation  and 
metamodel 

+n  :  number  of  observations 

++R3^^  :  the  coefficient  of  determination 
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TABLE  9. 


DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
HOLDOUT  GROUP  (C^^) 


Variable 

Mean 

S.D.* 

Min 

Max 

^Ai  o 

0.98 

0.03 

0.90 

1 . 02 

C  .  2 

0 . 98 

0.04 

0.83 

1 . 02 

r  .  2 

0 . 98 

0 . 04 

0.79 

1 . 02 

0.007 

0.007 

0.00022 

0 . 03 

£%*** 

0 . 70 

0.74 

0 . 022 

3.26 

n+  =  41, 

R,^^"  = 

0.8906 

*S.D. 

**£ 

***£% 


+n 

++R2^ 


standard  deviation 

absolute  error  of  Cm"  between  simulation  and 
metamodel 

percentage  error  of  Cm"  between  simulation  and 
metamodel 

number  of  observations 

the  coefficient  of  determination 


Departure  Processes 

The  departure -processes  module  estimates  the  mean  and 
squared  coefficient  of  variation  of  interdeparture  times. 
Based  on  the  flow  conservation  law,  if  the  V/C  ratio  is  less 
than  1,  the  mean  outflow  rate  should  be  equal  to  the  mean 
inflow  rate.  Therefore,  the  average  directional 
interdeparture  time  can  be  determined  from  the  corresponding 
interarrival  time: 


(19) 


The  coefficient  of  variation  of  interdeparture  times  is 
estimated  in  two  steps. 

Step  1 

The  coefficient  of  variation  of  interdeparture  times  is 
first  estimated  with  the  two  directions  combined.  Departure 
processes  with  generally  distributed  arrivals  and  service 
times  are  analyzed  by  using  Laplace  transforms.  The  use  of 
Laplace  transforms  for  derivations  in  queuing  theory  (which 
is  quite  frequent)  is  presented  in  texts  such  as  Kleinrock 
[21] .  Some  analytic  relations  obtained  in  this  dissertation 
are  shown  below  using  the  following  notation: 
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Let : 


tj,  o/ 

P 


mean  and  variance  of 
mean  and  variance  of 
mean  and  variance  of 
mean  and  variance  of 
V/C  ratio 


interarrival  times 
lock  service  times 
interdeparture  times 
lock  idle  times 


C  ^ 

"-A  ' 


C  ^ 
'-s  ' 


C  2 

'-D 


squared  coefficients  of  variation 
for  interarrival  times,  service 
times,  and  interdeparture  times 


fA(t)  ,  fs(t)  ,  f^{t)  ,  fi(t) 


probability  density  functions 
(pdf)  for  interarrival  times, 
lock  service  times, 
interdeparture  times,  and  lock 
idle  times 


F;(z)  ,  Fs*(z)  ,  Fo  (z)  ,  Fi*(z)  = 


Laplace  transforms  for 

f;,(t)  ,  fs(t)  ,  fD{t)  ,  fi(t) 


For  example,  for  interarrival  times,  the  Laplace  transform 
is  expressed  as 


Fa*  (2) 


=  I  fft(t)e-^^  dt 
0 


(20) 


The  departure  process  in  a  queuing  station  may  be 
analyzed  for  two  different  conditions:  with  and  without  a 
queue.  The  interdeparture -time  distribution  would  be  equal 
to  the  service-time  distribution  while  there  are  queues 
waiting  for  service.  However,  the  interdeparture  time  would 
be  equal  to  the  sum  of  the  idle  time  and  the  service  time 
while  there  is  no  queue.  Therefore,  the  Laplace  transforms 
for  the  interdeparture  time  distributions  can  be  represented 
as  follows: 

Fd  i  z)  \uith  queue  ~  Fs 
Fd  (without  queue  ~  Fi  i  z)  F5  (z) 


(21) 

(22) 


The  probability  of  having  a  queue  is  given  by  the 
volume/capacity  ratio  p  [21] .  Then,  the  probability  of  not 
having  a  queue  is  (1-p) .  Therefore,  the  Laplace  transform 
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for  the  interarrival -time  distribution  can  be  represented  by 
Eq.  23: 

Fd  {■Z)-(1  p)  Fq  {  z)  j  uichouc  queue*  1  wi  th  queue  (23) 

=  (1-p)  Fj*  (z)  Fg*  (z)  +pFs  (z) 

The  mean  of  a  distribution  can  be  represented  by  the 
negative  value  of  the  first  derivative  of  its  Laplace 
transform  when  z  equals  0.  Therefore,  the  average 
interarrival  time,  service  time,  interdeparture  time,  and 
idle  time  can  be  represented  by  Eqs.  24a  to  24d: 


dF^  (z) 
dz 


(24a) 


dFg  { z) 
dz 


(24b) 


3Fd  { z) 
dz 


(24c) 


~  ~ 


dFj  (z) 


(24d) 


The  variance  can  be  expressed  as  the  difference  between 
the  second  derivative  of  the  Laplace  transform  and  the 
squared  mean.  Eqs.  25a-25d  express  such  relations  for 
interarrival  time,  service  time,  interdeparture  time,  and 
idle  time  distributions: 


^Fft  (z) 


~  t  ^ 

z=0  ^A 


(25a) 


^Fo  (z) 


2  =  "  "S  I  _  -p  2 

c  _  _  -y-n  t'  a 


(25b) 


^^Fn  (z) 


2  =  "  I  _  -p  2 

n  _  *  I  •»=  n  n 


(25c) 


d^FUz) 


\z-o  -  t/ 


(25d) 


When  z  equals  0,  the  Laplace  transform  is  equal  to  1, 
producing  the  following  relations  for  interarrival -time , 
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service-time,  interdeparture -time,  and  idle-time 
distributions : 


f;(o) =1 

(26a) 

Fs*(0)  =1 

(26b) 

Fo  (0)  =1 

(26c) 

F;{0) =1 

(26d) 

Combining  Eqs . 


23,  24b,  24c,  24d,  26b,  and  26d  yields 
dF  *  ( z) 

- - =  -(  (1-p)  (-tj-ts)  ) 


(27) 


Due  to  flow  conservation,  if  the  V/C  ratio  is  less  than  1, 
then  the  average  interdeparture  time  would  be  equal  to  the 
average  interarrival  time: 


(28) 

Therefore,  Eqs.  27  and  28  can  be  combined: 

=  (1-p)  (Tj+t^) +pt5  =  (29) 

Since  p  =  tg/t^^,  Eq.  29  yields 

T,  (30) 


Combining  Eqs . 
and  30  yields 


23,  24b,  24d,  25b,  25c,  25d,  26b,  26c, 


(z) 

dz^ 


{l-p)a/  +  o/+(t^-t,,)  ts 


26d, 

(51) 


Dividing  Eq.  31  by  t^^,  we  can  obtain  the  following  relation 
for  the  squared  coefficient  of  variation  Cq^  : 

o  ^ 

=  (l-p)  +p+C/p^-p^  (32) 


In  the  special  case  where  the  arrival  process  is 
Poisson  distributed  and  the  service  times  are  exponentially 
distributed,  then  due  to  the  memoryless  property  of  the 


Poisson  distribution,  the  variance  of  idle  times  would  be 
equal  to  the  variance  of  interarrival  times.  Since  the 
interarrival  times  for  a  Poisson  process  are  exponentially 
distributed  and  since  the  mean  and  standard  deviation  of  an 
exponential  distribution  are  equal,  we  can  state  the 
following : 


o 


2 

I 


(33) 


a 


2 


S 


(34) 


Therefore,  in  this  special  case,  Eq.  31  can  be  simplified  to 


a 


2 

D 


(35) 


which  is  consistent  with  Burke's  theorem  [11] .  In  that 
theorem  Burke  proved  that,  when  the  arrivals  are  Poisson 
distributed  and  the  service  times  are  exponentially 
distributed,  then  the  departures  must  be  Poisson  distributed 
with  the  same  mean  and  variance  as  the  arrivals. 

The  main  difficulty  in  estimating  the  squared 
coefficient  of  variation  of  interdeparture  times  (Eq.  32) 
when  arrivals  and  service  times  are  generally  distributed  is 
determining  the  variance  of  the  lock  idle  times.  These 
depend  on  the  way  in  which  the  previous  busy  period 
terminated.  This  problem  may  be  bypassed  by  developing  a 
metamodel  for  directly  estimating  the  squared  coefficient  of 
variation  of  interdeparture  times. 

Before  developing  the  metamodel  for  estimating  the 
coefficient  of  variation  of  overall  interdeparture- time 
distribution,  the  stratum-specific  split-sample  scheme  was 
also  employed  to  split  the  database  of  Experiment  1  into  the 
two  groups  (training  and  holdout) .  For  developing  the 
metamodel,  there  were  120  observations  available  from  the 
database  of  Experiment  1.  Each  observation  was  obtained 
from  one  individual  lock  and  included  the  characteristics  of 
overall  interdeparture-time  distribution  and  corresponding 
overall  interarrival -time,  service-time  distributions,  and 
V/C  ratio  at  one  lock.  The  split  ratio  between  the  training 
and  holdout  groups  is  approximately  2:1.  Therefore,  the 
computer  randomly  assigned  79  observations  to  the  training 
group  and  41  observations  to  the  holdout  group. 

Following  the  approach  outlined  in  Section  4.1,  the 
following  metamodel  is  developed: 

=  0.9984  n  =  79  Sg  =  0.0058  =  0.83116 
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(36) 


0 . 207  +0 .795  (C/ (1-p)  -^p)  +1 . 001  (C/p2-p2) 

(0.065)  (0.066)  (0.0046) 

=  0.9984  n  =  79  =  0.0058  Py  =  0.83116 

The  metamodel  of  Eq.  36  was  developed  based  on  the 
structural  form  of  Eq.  32.  Eq.  36  was  originally  developed 
by  using  four  separate  variables:  (1)  C;^^(l-p),  (2)  p,  (3) 

Cg^p^,  and  (4)  p^.  However  very  high  correlations  were 
observed  between  Variables  1  and  2  and  Variables  3  and  4. 

The  correlation  coefficient  between  Variables  1  and  2  is  - 
0.9993  and  the  correlation  coefficient  between  Variables  3 
and  4  is  0.9222.  Moreove  the  coefficients  of  Variables  1 
and  2  were  almost  equal .  Variables  3  and  4  also  had  nearly 
equal  coefficients.  To  avoid  mult icollinearity  problems, 
Variables  1  and  2  were  combined  into  a  single  variable  while 
Variables  3  and  4  were  combined  into  a  second  variable  to 
develop  Eq.  36.  The  correlation  coefficient  between  the 
combined  variables  (1  and  2)  and  (3  and  4)  is  -0.1775,  which 
indicates  the  new  combined  variables  are  not  highly 
correlated.  It  is  noteworthy  that  the  sum  of  the  intercept 
(0.207)  and  the  parameter  of  Cp^^{l-p)+p  (0.795)  is 
approximately  equal  to  1 .  In  addition,  the  dependent 
variable  has  a  standard  error  of  0.0058,  which  is  only  0.7% 
of  its  mean. 

The  holdout  group  data  was  used  to  validate  the 
reliability  of  this  metamodel  and  the  cross-validation  test 
was  again  conducted.  The  descriptive  statistics  and  test 
results  for  the  training  group  and  the  holdout  group  are 
listed  in  Tables  10  and  11. 

This  metamodel  is  very  reliable  since  the  shrinjcage  is 
0.0029.  It  performs  especially  well  in  the  cross-validation 
test,  where  its  absolute  error  ranges  from  0.00014  to  0,03, 
with  a  mean  of  0.007,  and  its  absolute  percentage  error 
ranges  from  0.02  to  3.34,  with  a  mean  of  0.84. 

Since  the  mean  and  standard  deviation  of  an  exponential 
distribution  must  be  equal,  its  coefficient  of  variation 
must  be  1.0.  Thus,  for  the  special  case  of  an  M/M/l  queue: 

C/  =  1  (37) 


C/  =  1  (38) 
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TABLE  10.  DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
TRAINING  GROUP  {C^^) 


Variable 

!  Mean 

S.] 

D.* 

Min 

Max 

c  ^ 

0  . 

83 

0 

.  15 

0 

.39 

1 

.  02 

P  ^ 

0  . 

99 

0 

.  02 

0 

.  89 

1 

.  02 

C  2 

'“S 

0. 

53 

0 

.14 

0 

.34 

0 

.  81 

p  ^ 

0  . 

54 

0 

.27 

0 

.04 

0 

.  97 

£** 

0  . 

004 

0 

.  004 

0 

.  00009 

0 

.  02 

£%*** 

0. 

54 

0 

.46 

0 

.  02 

2. 

.30 

n*  =  79, 

Rt2+-' 

= 

0 . 9984 

*S.D. 

**£ 

***£% 


standard  deviation 

absolute  error  of  C^^  between  simulation  and 
metamodel 

percentage  error  of  Cj^^  between  simulation  and 
metamodel 

number  of  observations 

the  coefficient  of  determination 


TABLE  11.  DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
HOLDOUT  GROUP  (Cd^) 


Variable  Mean  S.D.*  Min 


*S.D.  :  standard  deviation 

**£  :  absolute  error  of  Cp^  between  simulation  and 

metamodel 

***£%  :  percentage  error  of  Cp^  between  simulation  and 

metamodel 

+n  :  number  of  observations 

++R,^  :  the  coefficient  of  determination 
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Substituting  Eqs .  37  and  38  into  Eq.  36,  the  latter  may  be 
simplified  as  follows: 

(39) 

=  0 . 207 +0 .795  (1-p+p)  +  (p2-p2)  =  0.207  +0.795  =  1.002  *  1.0 

This  result  is  also  consistent  with  Burke's  Theorem  [11]. 

The  parameters  and  structural  form  of  this  metamodel 
are  similar  to  those  of  Eq.  32,  which  was  analytically 
derived  for  G/G/l  queues.  In  addition,  its  standard  error 
is  extremely  tight  (Sg//iy  =  0.007)  and  it  satisfies  Burke's 
Theorem  very  closely  when  applied  to  the  special  M/M/1  case. 
It  may  be  concluded  that  the  good  performance  of  the 
metamodel  is  due  to  its  development  approach.  Its 
structural  form  was  based  on  queuing  theory  while  its 
coefficients  were  estimated  statistically  from  simulation 
results.  Such  a  metamodel  should  be  very  useful  for 
predicting  interdeparture  time  distributions  from  G/G/l 
queues  embedded  in  larger  systems,  such  as  series  and 
networks . 

Step  2 

The  coefficients  of  variation  of  directional 
interdeparture -times  for  upstream  and  downstream  traffic 
must  be  estimated.  For  this  purpose  the  following  metamodel 
was  developed: 

=  0. 518+0. 491C,,./  di  (40) 

(0.0056)  (0.0068) 

=  0.9710  n  =  158  Sg  =  0.013  /iy  =  0.9164 

where 

squared  coefficient  of  variation  of  directional 
interdeparture  times  for  Direction  j  and  Lock  i 

squared  coefficient  of  variation  of  directional 
interarrival  times  for  Direction  j  and  Lock  i 

:  squared  coefficient  of  variation  of 

interdeparture  times  for  Direction  j  and  Lock  i 

Before  developing  the  above  metamodel  (Eq.  40) ,  the 
stratum-specific  split-sample  scheme  was  also  employed  to 
split  the  database  of  Experiment  1  into  the  two  groups 
(training  and  holdout) .  For  developing  this  metamodel,  240 
observations  were  available  from  the  database  of  Experiment 
1.  Two  observations  (upstream  and  downstream  directional 
interdeparture-time  distributions)  could  be  obtained  from 
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one  individual  lock.  Each  observation  included  the 
directional  interdeparture-time  distribution  and 
corresponding  V/C  ratio,  the  overall  interdeparture- time, 
and  the  directional  interarrival-time  and  service-time 
distributions  at  one  lock.  The  split  ratio  between  the 
training  and  holdout  groups  was  approximately  2:1. 
Therefore,  the  computer  randomly  assigned  158  observations 
to  the  training  group  and  82  observations  to  the  holdout 
group . 

The  development  of  this  particular  metamodel  (Eq.  40) 
was  based  on  empirical  analysis.  The  Pearson  correlation 
coefficients  were  computed  and  the  dependent  variable 
(squared  coefficient  of  variation  of  directional 
interdeparture -time  distribution)  was  plotted  versus  each 
possible  influential  factor.  The  possible  influential 
factors  included  the  following: 


1. 

(squared  coefficient  of  variation 
interdeparture -time  distribution) , 

of 

overall 

2. 

(squared  coefficient  of  variation  of 
directional  interarrival-time  distribution) , 

3  . 

(squared  coefficient  of  variation 
interarrival-time  distribution) , 

of 

overall 

4  . 

Cg^  (squared  coefficient  of  variation 
times) , 

of 

service 

5  . 

V/C  ratio. 

6  . 

Ln(CD2)  , 

7. 

"“a  *“0  ' 

8  . 

(Cd2)0-5, 

9  . 

(Cg2)2^ 

10  . 

(C^2)0.5^ 

11 . 

Ca^Ln(CD^),  and 

12  . 

CD2Ln(Ca2) 

The  seventh  factor  in  the  above  list  was  selected 

because  it  yielded  the  highest  and  the  smallest  Sg . 

The  holdout  group  data  was  used  to  validate  the 
reliability  of  this  metamodel.  The  cross-validation  test 
was  again  performed.  The  descriptive  statistics  and  test 
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results  for  the  training  group  and  the  holdout  group  are 
shown  in  Tables  12  and  13 . 


The  shrinkage  (0.0007)  indicates  that  the  reliability 
of  this  metamodel  is  high.  In  addition,  the  results  of  the 
cross-validation  test  further  support  its  reliability  since 
the  absolute  error  ranges  from  0.00004  to  0.05,  with  a  mean 
of  0.009,  and  the  absolute  percentage  error  ranges  from 
0.005  to  5.41,  with  a  mean  of  0.96.  However,  the  maximum 
absolute  percentage  error  is  7%  in  the  training  group  and 
5.41  in  the  holdout  group.  Therefore,  this  metamodel  does 
not  always  predict  precisely  the  coefficient  of  variation  of 
directional  interdeparture-tima  distribution  and  could  stand 
further  improvement . 

Delay  Function 

The  delay  function  is  intended  to  estimate  the  average 
waiting  time  at  a  lock.  Marshall  [26]  has  tried  to  express 
the  variance  of  interdeparture  times  in  terms  of  the  average 
waiting  time.  In  this  study,  the  average  waiting  time  is 
expressed  in  terms  of  the  variance  of  interdeparture  times 
by  applying  Marshall's  formula.  An  exact  solution  for  the 
average  waiting  time  is  obtained  as  follows: 


2Tl{l-p) 


where 

W  :  the  average  waiting  time 

:  variance  of  interarrival  times 

:  variance  of  service  times 

:  variance  of  interdeparture  times 

t^  :  average  interarrival  time 
p  :  volume  to  capacity  ratio 


(41) 
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TABLE  12.  DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
TRAINING  GROUP  (C^^j 


Variable 


;dji- 

:ei 


e% 


n" 


=  158, 


Mean  S.D.*  Min  Max 


0  . 

.  92 

0 

.  08 

0 

.69 

1 . 

02 

0  . 

,  98 

0 

.  04 

0 

.  74 

1 . 

02 

0  , 

.  83 

0 

.15 

0 

.39 

1 . 

02 

0  , 

.  01 

0 

.  009 

0 

.00005 

0  . 

05 

1, 

.  08 

1 

.  08 

0 

.006 

7  . 

00 

=  0.9710 


*S.D. 

**£ 


***£% 

+n 

++R^^ 


standard  deviation 

absolute  error  of  between  simulation  and 

metamodel 

percentage  error  of  ^dji  between  simulation  and 
metamodel 

number  of  observations 

the  coefficient  of  determination 


TABLE  13 .  DESCRIPTIVE  STATISTICS  AND  TEST  RESULTS  FOR  THE 
HOLDOUT  GROUP  (C^^) 


Variable 

Mean 

S.D.* 

Min 

Max 

p  2 

^dDi2 

<^aji 

0.92 

0,98 

0.08 

0.03 

0.70 

0.79 

1.02 

1 . 02 

0.83 

0.15 

0.41 

1 . 02 

0.009 

0.0096 

0.00004 

0 . 05 

£%*** 

0 . 96 

1.03 

0.005 

5.41 

n-"  =  82, 

p  2  +  +  _ 

^2 

0.9703 

*S.D. 

**£ 

***£% 

+n 


standard  deviation 

absolute  error  of  C^j^^  between  simulation  and 
metamodel 

percentage  error  of  between  simulation  and 

metamodel 

number  of  observations 

the  coefficient  of  determination 
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For  the  general  queuing  system  the  difficulty  with  Eq. 
41  is  that  only  the  variances  of  interarrival  and  service 
times  are  known,  while  the  variance  of  interdeparture  times 
is  unknown.  This  difficulty  can  now  be  overcome  by  using 
the  metamodel  developed  in  this  study  to  estimate  the 
coefficient  of  variation  of  interdeparture  times  (Eq.  36). 

In  this  delay  function,  the  average  waiting  time 
increases  as  the  variance  of  interarrival  and  service  times 
increase  and  decreases  as  the  variance  of  interdeparture 
times  increases.  The  average  waiting  time  approaches 
infinity  as  the  volume/capacity  ratio  p  approaches  1.0. 

Algorithm  for  Two-Way  Traffic  Systems 

Delays  depend  on  the  interarrival,  interdeparture,  and 
service- time  distributions.  Therefore,  to  estimate  delays, 
we  need  to  know  in  advance  the  means  and  variances  of  the 
interarrival -time,  interdeparture-time,  and  service-time 
distributions.  For  two-way  traffic  systems  with  series  of 
G/G/1  queues  and  bi-directional  servers,  a  difficulty  arises 
in  identifying  the  variances  of  interarrival  and 
interdeparture  times.  The  variance  of  interarrival  times  at 
a  certain  lock  labeled  k  is  affected  by  the  interdeparture - 
time  distributions  from  adjacent  upstream  and  downstream 
locks.  The  interdeparture -time  distributions  at  adjacent 
upstream  and  downstream  locks  depend  on  their  interarrival - 
time  distributions,  which  are  affected  by  the 
interdeparture- time  distributions  from  Lock  k.  Hence,  the 
variances  of  interarrival  times  at  adjacent  locks  depend 
upon  each  other.  Therefore,  the  variances  of  interarrival 
times  cannot  be  determined  from  a  single  one-directional 
scan  along  a  series  of  queues.  For  example,  if  we  tried  to 
estimate  delays  by  scanning  from  upstream  toward  downstream, 
we  could  determine  at  Lock  k  the  variance  of  interarrival 
times  from  upstream,  but  the  variance  of  interarrival  times 
from  downstream  would  be  unknown  and  would  be  affected  by 
the  interdeparture-time  distribution  from  this  Lock  k.  To 
overcome  such  complex  interdependence,  an  iterative 
algorithm  is  proposed.  It  starts  scanning  in  one  direction 
while  using  some  initialized  assumed  values  for  the 
variances  of  interdeparture  times  from  the  opposite 
direction.  It  can  thus  sequentially  estimate  the 
interarrival  and  interdeparture- time  distributions  for  each 
lock.  Then,  the  scanning  direction  is  reversed  and  the 
process  is  repeated,  using  the  interdeparture-time 
distributions  for  the  opposite  direction  estimated  in  the 
previous  scan.  Alternating  directions,  the  scanning  process 
continues  until  the  specified  convergence  criteria  (usually 
the  squared  coefficient  of  variation  or  the  variances  of 
interdeparture  times)  computed  in  successive  iterations 


62 


converge.  Then  the  algorithm  stops  reestimating  the  arrival 
distributions  and  proceeds  to  estimate  delays. 

The  following  algorithm  is  designed  to  apply  the 
metamodels  in  estimating  single -chamber  lock  delays.  The 
notation  used  in  this  algorithm  is  defined  as  follows: 


average  interarrival  time  for  Direction  j  and 
Lock  i 


average  interdeparture  time  for  Direction  j  and 
Lock  k 


variance  of  interarrival  times  for  Direction  j , 
Lock  i  and  Iteration  n 


o 


2 

djkn 


variance  of  interdeparture  times  for  Direction  j, 
Lock  k  and  Iteration  n 


o 


2 

dlO 


variance  of  interdeparture  times  from  origin  node 


a 


2 

d2L*l 


variance  of  interdeparture  times  from  destination 
node 


V>vik 


^vik 

i 

j 

n 


Pi 

^i 


distance  between  Locks  i  and  k 

average  tow  speed  between  Locks  i  and  k 

standard  deviation  of  tow  speeds  between  Locks  i 
and  k 

lock  index 

direction  index  (1  =  downstream,  2  =  upstream) 
iteration 

average  interarrival  time  at  Lock  i 

variance  of  interarrival  times  at  Lock  i  and 
Iteration  n 

average  interdeparture  time  at  Lock  i 

variance  of  interdeparture  times  at  Lock  i  and 
Iteration  n 

variance  of  lock  service  times  at  Lock  i 
volume/capacity  (=V/C)  ratio  at  Lock  i 
average  inflow  rate  at  Lock  i 
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coefficient  of  variation  of  directional 
interarrival  times  for  Lock  i,  Direction  j,  and 
Iteration  n 


coefficient  of  variation  of  directional 
interdeparture  times  for  Lock  i,  Direction  j,  and 
Iteration  n 


'Si 


coefficient  of  variation  of  service  time  at  Lock 
i 


coefficients  of  variation  for 
interarrival  times  and  interdeparture 
times,  respectively,  at  Lock  i  and 
Iteration  n 


c-  :  constant  (assumed  value) 

£  :  total  number  of  Locks 


:  average  waiting  time  at  Lock  i 

m  :  indicator  of  scanning  direction  (1=  from  origin 

to  destination,  2=  from  destination  to  origin) 


The  required  inputs  of  this  algorithm  include  the  means 
and  coefficients  of  variation  for  service-time  distributions 
and  inflow  distributions,  distances  between  locks,  and  the 
means  and  standard  deviations  of  speed  distributions.  This 
algorithm  estimates  delays  as  well  as  means  and  coefficients 
of  varia'im  for  interarrival  and  interdeparture- time 
distribu  .  ons .  The  algorithm  consists  of  the  following 
steps : 


1. 


Compute  the  average  directional  interarrival  and 
interdeparture  times  for  each  lock: 

t  =t  =t  {ic=i-l,  if  7=1 

‘'dji  ^aji  ^djk  *ic=i+l,  if  7=2 


(15a) 


Compute  the  average  interarrival  time  for  each  lock: 


^ali *  ^a2i 


"Ai  -r 


^ali'^  ^a2i 


i=l,  .  .  .M 


(17a) 


3.  Estimate  the  coefficients  of  variation  for  interarrival 
and  interdeparture  times  at  each  lock. 

3.1  Set  n=l,  m=l 

3.2  Assume  initial  values  for  the  standard  deviations  of 
interdeparture  times  in  Direction  2,  at  Locks  2  through 
L-1 
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^  d2iO  ^  i  •  ^  ^  ^ 

3.3  Starting  from  Lock  1,  let  i  =  1 

3.4  Compute  the  standard  deviations  of  directional 
interarrival  times  at  Lock  i  using  the  metamodel 
expressed  in  Eq.  16: 


3.4.1  if  j=l  and  m=l 


Oai.n  =  Odi.-ln"0.0251  1  n  ( 1  ^  ) 

^vii-1 


(16a) 


3.4.2  if  j=2  and  m=l 


vi  i  *  1 


(16b) 


3.4.3  if  j=l  and  m=2 


<»aii.  =  Odii-ln-1^0-0251  Ind^^^iiii^^) 


(16c) 


3.4.4  if  j=2  and  m=2 

f*  vi  i  ♦  1 

3.5  Compute  the  squared  coefficients  of  variation  of 
directional  interarrival  times  at  Lock  i: 

o  ■  •  ^ 

r  2  _  /  '^ajin  \ 

'-ajin  '  —  > 

^aji 

3.6  Compute  the  squared  coefficient  of  variation  of 
combined  interarrival  times  at  Lock  i 


(16d) 


C  •  ^  =  0  179+0  41  (C  .  2  +  C 


(18a) 


3.7  Estimate  the  squared  coefficient  of  variation  for 
interdeparture  times  at  Lock  i  using  the  metamodel 
expressed  in  Eq.  36: 

Coin  =  0. 207+0. 795  (C^,„dl-p,)+p,)+1.001(C5/p^p^)  (36a) 

3.8  Compute  the  variance  of  interarrival  times  at  Lock  i 


T  .  ^  =  P  2  ^ 

Ain  ^Ain  ^Ai 
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3.9  Compute  the  variance  of  interdeparture  times  at  Lock  i 


o 


2 

Din 


c  ^  t  ^ 

'~Din  ^Di 


(45) 


3.10  Compute  the  squared  coefficients  of  variation  of 

directional  interdeparture  times  at  Lock  i,  using  the 
relation  developed  in  Eq.  40: 

-  0. 518+0. 491C,.,„2  c^.^2  (4Ca) 


3.11  Compute  the  standard  deviations  of  directional 
interdeparture  times  at  Lock  i: 

^djin  ~  ^djin  ^dji 


(46) 


3.12  Repeat  Steps  3.4  -  3.11  for  i  =  2,  ...L 

3.13  Set  n  =  n+1,  m=2 

3.14  Starting  from  Lock  L,  let  i  =  L,  and  repeat  Steps 

3.4  -  3.12  for  i  =  L,  (L-1) ,  ...1 

3.15  Set  n  =  n+1,  m=l 

3.16  Repeat  Steps  3.4  -  3.14 

3.17  If  the  following  condition  is  satisfied,  then  go  to 
Step  4.  Otherwise,  go  to  Step  3.15 


If  2  _  f  2  1 

I  ^djin _ ^djin-1  I 


C  ^ 


^  0.001 


i-1,  ...L,  j-1,2 


(47) 


4 .  Estimate  the  average  waiting  times  using  the  formula 
expressed  in  Eq.  41: 


= 


'Ain 


+2o 


S-l 


'Din 


2t„,(l-p,) 


(41a) 


Algorithm  for  One-Way  Traffic  Systems 

Although  the  numerical  method  was  originally 
developed  for  two-way  traffic  systems,  with  a  few 
simplifications  this  methed  can  be  adapted  for  the  more 
generally  encountered  systems  with  one-directional  servers. 
The  one-directional  systems  may  be  treated  as  a  special  case 
of  the  two-directional  systems.  The  one -direct ional 
algorithm  should  perform  better  since  the  interarrival  time 
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distributions  will  be  affected  by  the  interdeparture  time 
distributions  from  upstream  only  and  are  not  subject  to 
circular  interdependence.  Therefore,  the  interarrival - t ime 
distributions  may  be  determined  in  a  single  one -direct ional 
scan  without  any  iteration.  The  notation  used  in  this 
algorithm  for  one-way  traffic  is  defined  as  follows: 

D  ,.  distance  between  Locks  i  and  k 

ik 


vik 


average  tow  speed  between  Locks  i  and  k 


^vik 


standard  deviation  of  tow  speeds  between  Locks  i 
and  k 


Pi 

^i 


the  average  interarrival  time  at  Lock  i 
variance  of  interarrival  times  at  Lock  i 
average  interdeparture  time  at  Lock  i 
variance  of  interdeparture  times  at  Lock  i 
variance  of  interdeparture  times  from  origin  node 
variance  of  lock  service  times  at  Lock  i 
V/C  ratio  at  Lock  i 
average  inflow  rate  at  Lock  i 


2  2  2 
^Ai  •  ^Di  •  ^Si 


squared  coefficients  of  variation  for 
interarrival  times,  interdeparture 
times,  and  service  times  respectively, 
at  Lock  i . 


:  total  number  of  locks 

Vl^  :  the  average  waiting  time  at  Lock  i 

This  algorithm  has  the  same  inputs  and  outputs  as  the 
algorithm  for  two-way  systems.  The  algorithm  consists  of 
the  following  steps: 

1.  Compute  the  average  interarrival  time  and 
interdeparture  time  for  each  lock: 

^Di  ~  ^Ai  ~  ^Di-1  -^“1'  •  ■  -  L  (15b) 
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2.  Estimate  the  squared  coefficients  of  variation  and, 
from  them,  the  variances  of  interarrival  and 
interdeparture  time  distributions  for  each  lock. 

2.1  Starting  from  Lock  1,  let  i=l 

2.2  Compute  the  standard  deviation  of  interarrival  times 
at  Lock  i  using  the  metamodel  developed  for  two-way 
traffic  (Eq.  16 ) : 


.1  +  0 . 0  2  5 1 1  n  ( 1  +  ) 


(16e) 


2.3  Compute  the  squared  coefficient  of  variation  for 
interarrival  times  at  Lock  i 


o  ^ 

^Ai 


(48) 


2.4  Estimate  the  squared  coefficient  of  variation  for 
interdeparture  times  at  Lock  i  using  the  metamodel 
developed  for  two-way  traffic  (Eq.  36)  : 

=  0 .207+0 .795  (1-p^) +p^) +1. 001  (C^/pi-pi)  (36b) 


2.5  Compute  the  variance  of  interdeparture  times  at  Lock  i 


'Di 


2  =  (7  .2(-  2 
^Dl  '-Dl 


(49) 


2 . 6 

3  . 


Repeat  steps  2.2  -  2.5  for  i  =  2, 


.L 


Estimate  average  waiting  times  using  the  relation 
developed  for  two-way  traffic  (Eq.  41) : 


W,  = 


^Ai^°Si  -^Di 

2t^,{l-p,) 


(41b) 
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CHAPTER  5 

VALIDATION  OF  NUMERICAL  METHOD 


Comparison  of  Nximerlcal  and  Simulated  Results 

An  experiment  is  conducted  to  test  how  well  the 
numerical  method  duplicates  the  results  of  simulations. 

Eight  three-lock,  two-directional  systems  are  tested.  The 
controlled  variables  in  this  experiment  include  the  V/C 
ratio  p,  the  variance  of  lock  service  times,  inflow  rate, 
distance  between  locks,  and  tow  speed.  Values  for  these 
variables  are  generated  randomly  and  uniformly  between 
specified  bounds  which  cover  the  ranges  of  values  where 
waterway  delays  can  be  significant.  In  this  test  the  ranges 
of  the  controlled  variables  are  0.01  to  0.89  for  the  V/C 
ratio,  0.0007  to  0.3332  for  the  variance  of  lock  service 
times,  6.0  to  28.5  tows  per  day  for  the  inflow  rate,  5  to  60 
miles  for  the  distance  between  locks,  108  to  325  miles  per 
day  for  the  average  tow  speed,  and  33.84  to  101.52  miles  per 
day  for  the  standard  deviation  of  tow  speeds.  Table  14 
lists  the  values  of  controlled  variables  for  each  system  in 
this  experiment.  For  this  system  the  variances  of  the 
interdeparture  times  converged  within  0.1%  for  every  lock 
and  direction  in  no  more  than  6  iterations.  The  results  are 
shown  in  Tables  15  and  16. 

The  performance  of  the  numerical  method  could  be 
discussed  in  terms  of  the  scanning  process  and  the  average 
waiting  time.  The  performance  of  scanning  processes  is  the 
overall  performance  of  arrival,  departure  and  iterative 
scanning  processes.  The  performance  indicators  are  and 
£-q,  which  directly  measure  the  variance  deviations  of 
interarrival  times  and  interdeparture  times,  respectively, 
between  numerical  and  simulated  results.  The  average 
absolute  6^  is  3.57%  for  these  8  systems.  For  17  of  24  locks 
in  these  8  three- lock  systems  is  less  than  4%.  The 
average  absolute  is  1.80%  for  these  8  systems.  For  23  of 
the  24  locks  £p  is  less  than  4%.  The  results  show  that  the 
scanning  process  closely  approximates  the  simulation 
estimates  for  the  variances  of  interarrival  and 
interdeparture  times. 

The  deviation  of  average  waiting  times  are  due  to  the 
scanning  process  and  the  delay  function.  To  clearly  define 
the  deviation  of  average  waiting  times,  the  following 
notations  are  used: 
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Wg^n,:  simulated  waiting  times 

Wnum=  waiting  times  estimated  with  numerical  method 

Wnutn2=  waiting  times  estimated  by  Eq.  41  with  the 
simulated  variances  of  interarrival  and 
interdeparture  times 

AW:  representing  the  total  deviation 

between  the  numerical  method  and  simulation 

AW2 :  representing  the  deviation  due  to  the 

delay  function 

AWl :  AW-AW2,  representing  the  deviation  due  to  the 
scanning  process 

AW/Wgj^^*100%,  representing  the  total  percentage 
deviation  between  numerical  and  simulated  results 

£^2’  AW2/Wgj^^*100% ,  representing  the  percentage 

deviation  due  to  the  delay  function 

^wi •  ^w"^W2'  representing  the  percentage  deviation  due 

to  the  scanning  process 

A  question  arises  about  the  deviation  between  and 

^num2  AW2)  .  The  obtained  by  applying  Eq.  41 

while  the  input  variances  of  interarrival  and  interdeparture 
times  are  simulated  results.  Since  the  inputs  are  simulated 
results  and  Eq.  41  is  an  exact  formula  for  average  waiting 
times,  why  is  there  some  deviation  between  Wg^^  and 
This  test  assumed  the  simulated  results  exactly  represent 
the  true  values.  However,  simulation  is  a  kind  of 
experiment  and  thus  its  results  are  always  subject  to 
certain  errors.  More  replications  and  longer  simulation 
periods  can  reduce,  but  not  completely  eliminate,  the 
errors.  That  is  why  there  is  some  deviation  between  Wg^^^ 
and  W^^^2  •  Table  15  shows  that  some  locks  have  quite  high 
£„2'  Table  16  shows  that  for  21  of  the  24  locks  absolute 

AW2  values  are  below  0.1  hr,  and  for  10  locks  the  absolute 
AW2  values  are  below  0.01  hr.  It  may  be  concluded  that  the 
delay  function  produces  results  consistent  with  simulated 
ones . 


AWl  and  indicate  the  deviation  and  percentage 
deviation  of  average  waiting  times  from  the  scanning 
process.  Table  15  shows  that  10  of  the  24  locks  have  £„^ 
values  above  10%.  However,  Table  16  shows  that  21  of  the  24 
locks  have  absolute  AWl  values  below  0.1  hr,  and  8  locks 
have  values  below  0.01  hr.  Therefore,  the  deviations 
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attributable  to  the  scanning  process  are  quite  minor  in  most 
cases . 

The  total  percentage  deviation  of  waiting  times  in 
Table  15  shows  that  8  of  the  24  locks  have  absolute  values 
above  10%.  The  total  deviation  of  waiting  times  AW  in  Table 
16  shows  that  21  of  the  24  locks  have  absolute  values  below 
0.1  hr,  and  11  locks  have  values  below  0.01  hr.  Therefore, 
the  numerical  method  works  satisfactorily.  In  general, 
these  results  indicate  that  the  numerical  method  may  be  used 
to  screen  alternatives  and  greatly  reduce  the  number  of  lock 
improvement  combinations  that  have  to  be  evaluated  by  the 
more  detailed  microscopic  simulation  model. 

Three  types  of  errors  arise  in  the  numerical  method: 
from  simulation,  from  the  metamodeling  procedure,  and  from 
the  iterative  scanning  process.  The  errors  from  simulation 
could  be  reduced  by  increasing  the  number  of  replications 
and  duration  of  simulated  periods.  The  errors  from  the 
metamodeling  procedure  could  be  reduced  by  collecting  more 
data  (that  is,  increasing  the  data  base) .  Such  improvements 
could  increase  the  accuracy  of  the  numerical  method. 

To  apply  the  numerical  method  to  more  general  networks 
of  queues,  it  is  desirable  to  build  a  new  data  base  for 
redeveloping  or  validating  metamodels  since  the  applicable 
ranges  of  the  current  one  are  mainly  appropriate  for 
waterways . 
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TABLE  14.  VALUES  OF  CONTROLLED  VARIABLES 


System  Lock  V/C  Dist  Speed 

tows/day  mi  mi/day 

Mean  S.D.*^ 


1 

1 

6 . 0 

0 .01 

5 

270 

85 

0 . 0007 

2 

6.0 

0.07 

5 

270 

85 

0 . 0360 

3 

6 . 0 

0 . 17 

5 

270 

85 

0 . 1897 

2 

1 

12 . 0 

0 . 15 

5 

325 

102 

0 . 0309 

2 

12 . 0 

0 . 34 

5 

325 

102 

0 . 1620 

3 

12 . 0 

0.25 

5 

325 

102 

0 . 0915 

3 

1 

18 . 0 

0.22 

5 

108 

34 

0 . 0309 

2 

18 . 0 

0.03 

5 

108 

34 

0 . 0006 

3 

18 . 0 

0 . 50 

5 

108 

34 

0 . 1618 

4 

1 

24 . 0 

0 . 50 

5 

162 

51 

0 . 1883 

2 

24 . 0 

0.29 

5 

162 

51 

0 . 0646 

3 

24 . 0 

0.67 

5 

162 

51 

0.3330 

5 

1 

27.0 

0.75 

10 

108 

34 

0.2271 

2 

27 . 0 

0 . 57 

10 

108 

34 

0 . 1279 

3 

27 . 0 

0 .89 

10 

108 

34 

0.3167 

6 

1 

27.0 

0.75 

20 

216 

68 

0.1616 

2 

27.0 

0.57 

20 

216 

68 

0 .0909 

3 

27 . 0 

0.89 

20 

216 

68 

0 . 2259 

7 

1 

28 . 5 

0 . 60 

5 

325 

102 

0 . 1557 

2 

28 . 5 

0 . 05 

5 

325 

102 

0 . 0011 

3 

28 . 5 

0 . 80 

5 

325 

102 

0.2738 

8 

1 

28 . 5 

0 .35 

60 

162 

51 

0 . 0645 

2 

28 . 5 

0.60 

60 

162 

51 

0 . 1882 

3 

28 . 5 

0 . 80 

60 

162 

51 

0.3332 

*1  X  :  two-way  flow  rate 

*2  :  variance  of  service  times 

*3  S.D.:  standard  deviation 
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TABLE  15.  COMPARISON  OF  NUMERICAL  AND  SIMULATED  RESULTS  (1) 


Sys*^ 

Lock 

hr 

AW*^ 

hr 

£ 

o, 

o 

£ 

‘-wi 

o, 

"o 

£ 

^W2 

O, 

"O 

£ 

% 

£ 

% 

1 

1 

0 . 0003 

.  0015 

484 . 95 

584 . 95 

-100 . 00 

1 . 20 

1 . 09 

2 

0 . 0153 

.  0054 

35 . 57 

67.94 

-32 .37 

1 . 90 

1.42 

3 

0 . 0989 

.  0042 

4.22 

8 . 17 

-3 . 95 

0 . 95 

0 . 63 

2 

1 

0 . 0334 

- .0010 

-3.08 

47.43 

-50.51 

3 .38 

1 . 98 

2 

0.2316 

.  0047 

2 . 04 

14 . 87 

-12.83 

3 .81 

1 . 55 

3 

0 . 1139 

.0009 

0 .81 

25 .39 

-24 . 58 

3 . 20 

0 . 99 

3 

1 

0 . 0542 

.  0030 

5 . 53 

12.45 

-6 . 92 

3 .20 

2.45 

2 

0 . 0008 

- .0008 

-100.00 

0 . 00 

-100 . 00 

3 . 74 

3 . 94 

3 

0.4621 

.0155 

3.34 

4 . 63 

-1.29 

3 . 14 

1 . 82 

4 

1 

0.4355 

.  0155 

3 . 55 

6.40 

-2 . 85 

2 . 57 

-0 . 23 

2 

0 . 0962 

.  0094 

9.79 

28 . 73 

-18 . 94 

5.75 

1 .59 

3 

1.2028 

.  0016 

0.13 

4 . 18 

-4 . 05 

2 . 51 

-0 . 91 

5 

1 

1 .3926 

.1155 

8.30 

-0,41 

8 . 71 

-0 . 78 

-0 . 62 

2 

0.3901 

.  0413 

10.59 

-0 . 59 

11 . 18 

-0 . 15 

0 . 11 

3 

4 . 9837 

- .0868 

-1.74 

0.27 

-2 . 01 

-0 . 15 

-0 . 79 

6 

1 

1.2203 

.1339 

10 . 97 

-0 . 67 

11 . 64 

-0 . 50 

-0 . 07 

2 

0.3286 

.  0371 

11.29 

-4 . 61 

15.90 

-0.53 

1.36 

3 

4.4608 

- .0278 

-0.62 

1.09 

-1 . 71 

0.73 

-0.93 

7 

1 

0 . 5430 

.  0615 

11.32 

10.57 

0 . 75 

6 .59 

0 . 85 

i 

2 

0 . 0012 

-  .0012 

-100.00 

0.00 

-100 . 00 

10.73 

10 . 82 

3 

2 . 0874 

.  0391 

1.87 

5.42 

-3 . 55 

5 . 14 

-0 . 68 

8 

1 

0 . 1372 

.  0227 

16.57 

16 . 11 

0.46 

6 . 51 

3 . 23 

2 

0.6381 

.  0722 

11.32 

10.30 

1 . 02 

9 . 07 

3 . 03 

3 

2 . 3165 

.  1134 

4 . 89 

6 . 79 

-1 . 90 

9.46 

2 . 18 

*1  Sys  :  System 

*2  :  simulated  waiting  time 

*3  AW  :  where  is  waiting  times  estimated  with 

numerical  method 
*4  £„  :  AW/W3i^*100% 

*5  '  ^W~^W2 

*6  £„2  :  (Wjj^m2-Wsitn) /Wsim*lOO^.  where  is  waiting  times 

estimated  by  Eg.  41  with  simulated  variances  of 
interarrival  and  interdeparture  times 
*7  £a  =  where  is  variance  of 

interarrival  times  estimated  by  numerical  method; 
^A^aim  i®  simulated  variance  of  interarrival  times 
*8  £d  :  /(^D^sim*100%,  where  is  variance  of 

interarrival  times  estimated  by  numerical  method; 
ao^sitn  is  simulated  variance  of  interarrival  times 
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TABLE  16.  COMPARISON  OF  NUMERICAL  AND  SIMULATED  RESULTS  (2) 


Sys*^ 

Lock  Wgi^*2 

hr 

AW*3 

hr 

AWl*^ 

hr 

AW2*^ 

hr 

1 

1 

0 . 0003 

.  0015 

0 . 0018 

-0 . 0003 

2 

0 . 0153 

.  0054 

0 . 0104 

-0 . 0050 

3 

0 . 0989 

.  0042 

0 . 0031 

-0 . 0039 

2 

1 

0 . 0334 

- . 0010 

0 .0159 

-0.0169 

2 

0.2316 

.  0047 

0 . 0344 

-0 . 0297 

3 

0 . 1139 

.  0009 

0 . 0289 

-0.0280 

3 

1 

0.0542 

.0030 

0 . 0067 

-0.0037 

2 

0 . 0008 

- . 0008 

0 . 0000 

-0 . 0008 

3 

0.4621 

.  0155 

0 . 0214 

-0 . 0059 

4 

1 

0.4355 

.  0155 

0 . 0279 

-0.0124 

2 

0 . 0962 

.  0094 

0 . 0276 

-0 . 0182 

3 

1.2028 

.  0016 

0 . 0504 

-0.0488 

5 

1 

1 .3926 

.  1155 

-0 . 0058 

0.1213 

2 

0.3901 

.  0413 

-0.0023 

0.0436 

3 

4 . 9837 

- . 0868 

0.0135 

-0.1003 

6 

1 

1.2203 

.  1339 

-0.0081 

0 . 1420 

2 

0.3286 

.0371 

-0.0151 

0.0522 

3 

4.4608 

- . 0278 

0 . 0487 

-0.0765 

7 

1 

0 . 5430 

.  0615 

0 . 0574 

0.0041 

2 

0 . 0012 

- . 0012 

0.0000 

-0.0012 

3 

2 . 0874 

.  0391 

0.1132 

-0.0741 

8 

1 

0 . 1372 

.  0227 

0 . 0221 

0.0006 

2 

0 . 6381 

.  0722 

0.0657 

0.0065 

3 

2 .3165 

.  1134 

0 . 1573 

-0.0439 

*1  Sys 
*2  W  • 

'’sim 

*3  AW 

*4  AWl 
*5  AW2 


System 

simulated  waiting  time 

Wnum'Wgin,,  where  is  waiting  times  estimated 

with  numerical  method 

AW-AW2 

^num2"'^siin'  where  Wj„.^2  waiting  times  estimated 
by  Eq.  41  with  simulated  variances  of  interarrival 
and  interdeparture  tim.es 
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APPENDIX  A 


List  of  Variables 


Variable  Description 

ANB(i) :  Average  number  of  barges  for  0-D  pair  i,  i=l,NP 

BIAS(i) :  Chamber  assignment  policy  for  favoring  the  use  of 

main  chamber  at  lock  i,  i=l,NL,  (min) 

DIST(i):  Link  length,  i=l,NL,  (miles) 

DISTL(i) :  Distance  between  in-node  to  lock  on  link  i, 
i=l,NL,  (miles) 


DTEX(i) :  Mean  of  dwell  time  for  0-D  pair  i,  i=l,NP,  (days) 

DTSTD(i) ;  Standard  deviation  of  dwelT  time  for  0-D  pair  i, 
i=l,NP,  (days) 

EXS :  Mean  of  tow  speed,  (miles/day) 

FRB(i, j) :  The  lower  bound  of  interval  j  in  the  tow  size 
distribution  table  for  0  D  pair  i,  i=l,NP, 
j=l,NBI+l 

FRI(i) :  The  lower  bound  of  interval  i  for  trip  generation 

table,  i  =  l,NIl4-l 

FRS (i , j , k, 1 ) :  The  lower  bound  of  interval  1  for  lock 

service  time  table  for  lock  i,  chamber  j,  cut 
k,  i=l,NL,  j=l,NPL(i),  k=l,ML(i,j),  L=1,NSI+1 

FSTA(i,j) :  Average  stall  frequency  at  lock  i,  chamber  j, 

(#/yr) 

ID(i) :  Index  of  destination  for  0-D  pair  i,  i-l,NP 

IDBG:  =1,  to  actuate  the  debugging  function 

IDBGl :  =1,  to  debug  the  sequence  of  time  advance  and  event 
type 

IDBG2 :  =1,  to  debug  the  lock  event 
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IDBG3 :  =1,  to  debug  the  lock  event  only  at  certain  locks 

ILBG(i) :  =1,  to  actuate  the  debugging  function  at  lock  i, 

i=l,NL 

INS:  Specifying  the  starting  inventory  level  INS  while 

printing  detail  inventory  information 

IO(i) :  Index  of  origin  for  0-D  pair  i,  i=l,NP 

IGF:  =1,  to  plot  inventory  level 

lOI :  =1,  to  print  input  data 

lOM;  =1,  to  provide  the  inventory  information 

lOU:  =1,  to  plot  cumulative  delivery  and  consumption 

IP:  Indicator  for  printing  wait  time  data 

=0,  no  printing 

=1,  by  lock  and  direction  in  different  files 
=2,  in  one  file 

IXB(i,  j)  :  Random  r’-imber  seed  of  tow  size  (number  of  barges) 
for  0-D  pair  i,  direction  j,  i=l,NP,  j=l,2 

IXC(i) :  Random  number  seed  for  coal  consumption  at  node  i, 

i=l,NN 

IXL(i, j,k) :  Random  number  seed  for  lock  service  time  of  k 

cuts  at  lock  i,  chamber  j,  i=l,NL, 
j=l,NPL(i) ,  k=l,ML{i, j) 

IXS :  Random  number  seed  for  tow  speed 

IXSD(i,j) :  Random  number  seed  for  stall  duration  at  lock 

i,  chamber  j,  i=l,NL,  j=l,NPL(i) 

IXST(i, j) ;  Random  number  seed  for  stall  frequency  at 

lock  i,  chamber  j,  i=l,NL,  j=l,NPL(i) 

IXT(i,j):  Random  number  seed  for  inter-trip-generation- time 
for  0-D  pair  i,  direction  j,  i=l,NP,  j=l,2 

KSG :  Indicator  of  trip  generation  distribution 

=1,  empirical  distribution 
=2,  Poisson  distribution 
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KST:  Indicator  of  lock  service  time  distribution 

=1,  empirical  distribution 
=2,  uniform  distribution 
=3,  exponential  distribution 

ML{i, j) :  Types  of  cuts  at  lock  i,  chamber  j,  i=l,NL, 
j=l,NPL{i) 

NBI :  Total  number  of  intervals  in  the  distribution  table 

for  the  number  of  barges 

NCB(i, j,k) ;  Upper  limit  on  cut  size  k,  at  lock  i,  chamber 
j,  i=l,NL,  j=l,NPL(i),  k=l,ML(i,j) 

NC :  0-D  pairs  1-NC  are  coal  traffic 

NCTL(i) :  Lower  bound  of  number  of  coal  barges  per  coal  tow 

for  coal  0-D  pair  i,  i=l,NC 

NCTU(i) :  Upper  bound  of  number  of  coal  barges  per  coal  tow 

for  coal  0-D  pair  I,  1=1, NC 

NI(i) ;  Index  of  in  node  for  link  i,  i=l,NL 

Nil:  Total  number  of  intervals  in  the  trip  generation 

table 

NIS:  Number  of  different  starting  inventory  levels 

NL;  Total  number  of  locks/links 

NN :  Total  number  of  nodes 

NO(i) :  Index  of  out  node  for  link  i,  i=l,NL 

NP:  Total  number  of  0-D  pairs 

NPL(i) :  Total  number  of  chambers  in  lock  i,  i=l,NL 

NRAN :  Index  of  random  number  set 

NSI :  Total  number  of  intervals  in  lock  service  time 

distribution  table 

NTP :  Total  number  of  simulation  time  periods 
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NTR:  Index  of  trip  rate  set 

PLOAD:  Barge  payload,  (short -tons) 

RCON(i, j) :  Average  consumption  rate  for  simulation  time 

period  i,  at  node  j,  i=l,NTP,  j=l,NN,  (short- 
tons/day) 

RCT(i) :  Coal  barge  fraction  of  a  tow  for  0-D  pair  i, 

i=l,NC 

RCTL(i) :  Lower  bound  ratio  of  RCT(i),  i=l,NC 

RCTU(i) :  Upper  bound  ratio  of  RCT(i),  i=l,NC 

RL(i) ;  Lower  bound  ratio  of  average  consumption  rate  at 
node  i,  i=l,NN 

RR(i) :  Upper  bound  ratio  of  average  consumption  rate  at 

node  i,  i=l,NN 

RSTA(i,j) :  Stall  duration  at  lock  i,  chamber  j,  i=l,NL, 

j=l,NPL(i) 

RSTD(i,j,k):  Ratio  to  tighten  the  standard  deviation  for 

lock  service  time  distribution,  i=l,NL, 
j=l,NPL{i),  k=l,ML(i,j) 

RV;  Ratio  of  backhaul  speeds  to  linehaul  speeds 

SFST(i,j) :  jth  starting  inventory  level  at  node  i, 

i=l,NN,  j=l,NIS 

SRL(i) :  Lower  bound  ratio  of  TLO(i),  i=l,NL 

SRR(i) :  Upper  bound  ratio  of  TLO(i),  i=l,NL 

STDS:  Standard  deviation  of  tow  speeds 

TLO(i,j,k) :  Average  lock  service  time  for  cut  type  k,  at 

lock  i,  chamber  j,  i=l,NL,  j=l,NPL(i), 
k=l,ML(i,j)  (days/tow) 

TN(i, j) :  Trip  rate  for  0-D  pair  i,  simulation  time  period 
j,  i=l,NP,  j=l,NTP,  (tow  trips/day) 

TSTOP(i) :  End  of  simulation  time  period  i,  i=l,NTP,  (day) 
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TUT: 


Specified  time  interval  for  printing  inventory- 
levels,  (days) 


TWT:  Specific  time  period  to  provide  queue  length  data, 

(day) 

Input  Format 

Model  Input 

The  model  input  is  divided  into  13  sets  and  is 
described  below: 

1.  Basic  model  parameters:  (1  line) 
NN.NL.NP.NTP.NC.NIS. KST, KSG. IP  (1216) 

2.  Links/Loc)cs  characteristics:  (  (4  +  3*NPL  (  i)  )  *NL 
lines) 

(1)  NI (i) .NQ(i) ■NPL(i)  (1216) 

(2)  DIST(i) .DISTL(i)  (10F7.2) 

(3)  FSTA(i. i ) . i=l.NPL(i)  (10F7.2) 

(4)  RSTA(i.  i  )  .  i=l.NPL(i)  (10F7.'’) 

(5)  Chamber  characteristics:  (3*WPL(i)  lines) 

a.  ML  ( i  .  -i )  (1216) 

b.  TLO(i.  i  .k)  .]c=l.ML(i.  i  )  (8F8.5) 

c.  NCB(i.i.k)  .k=l.ML(i.-i)  (1216) 

j=l,NPL(i) 

i=l,NL 

3.  Traffic  demand:  (3*NP  lines) 

(1)  IO(i) . ID(i)  (1216) 

(2)  DTEX(i) .DTSTD(i)  (10F7.2) 

(3)  ANB (i)  (10F7.2) 

i=l,NP 

4.  Coal  tow  characteristics:  (3  lines,  if  NC>0) 

(1)  RCT(i) ■ i=l.NC  (10F7.2) 

(2)  NCTL(i) . i=l.NC  (1216) 

(3)  NCTU(i) .i=l.NC  (1216) 

5.  Speeds:  (2  lines) 

(1)  EXS . STDS  (10F7.2) 

(2)  RV  (10F7.2) 

6.  Simulation  time  periods  (NTP  lines) 

TSTOP  (i)  (10F7.2) 

i=l,NTP 


SS 


7  . 


8  . 


9  . 


10 


11 
12  , 


13 


Lock  service  time  distributions: 

(1)  if  KST=1  (NL*NPL(i) * (NSI+1) +1  lines) 

a.  (1216) 

b.  FRS(i.i .k.l) .k=l.ML(i.i)  (6F12.4) 

1=1,NSI+1 

j=l,NPL(i) 

i=l,NL 

(2)  if  KST=2  (NL  lines) 

a.  SRL(i) .SRR(i)  (10F7.2) 
i=l,NL 

Trip  generation  distributions: 

(1)  if  KSG=1  (NII+2) 

a.  Nil  (1216) 

b.  FRI (i)  (6F12.4) 

i=l,NII+l 


(NP* (NBI+1 ) +1  lines] 


size  distributions; 
(1)  NBI  (1216) 

(4)  FRB  ( i  .  -i  1  (6F12.4) 

j=l,NBI+l 
i=l,NP 


Inventories:  (2*NN+NTP+3) 

(1)  if  NC>0 

a.  SFST(i.i) .  i=l.NIS  (7F10.1) 
i=l,NN 

b.  PLOAD  (2F10.1) 

c.  RCON(i.  -i )  ■  i=l.NN  (10F7.1) 
i=l,NTP 

d.  RL ( i ) . RR ( i )  (10F7.2) 

i=l,NN 

e.  I^  (1216) 

(2)  TUT  (10F7.2) 

BIAS(i) .i=l.NL  (10F7.2)  (1  line) 

Debugging  function:  (2  or  3  lines) 

(1)  lOI. lOM. lOU. lOF  (1216) 

(2)  IDBG . IDBGl . IDBG2 . IDBG3 . IDBG4 . IDBG5 . IDBG6 
(1216) 

(3)  if  IDBG3=1 

a.  ILBG(i) . i=l.NL  (1216) 

if  KST=1: 

(1)  RSTD(i.i.k) .k=l.ML(i.i)  (6F12.4) 
j=l,NPL(i, j) 
i=l,NL 


86 


Random  Number  Seeds  Input 

1.  NRAN  (7111)  (1  line) 

2.  ( (NPL(i) +2) *NL  lines) 

(1)  IXSTd.  i)  .  i=l.NPL(i)  (7111) 

(2)  IXSD(i.i) .i=l.NPL(i)  (7111) 

(3)  IXL(i. i .k) .k=l.ML(i. i )  (7111) 

j=l,NPL(i) 

i=l,NL 

3.  (3*NP  lines) 

(1)  IXDT(i)  (7111) 

(2)  IXT(i.i) ,IXB(i.i)  (7111) 
j=l,2 

i=l,NP 

4.  I}^  (7111) 

5.  if  NC>0  (NN  lines) 

(1)  IXC(i)  (7111) 

6.  TWT  (10F7.2)  (1  line) 

Trip  Volume  Input 

1.  NTS  (1216) 

2.  TN(i.i) .i^l.NTP  (12F6.2) 
i=l,NP 
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APPENDIX  B 


Residual  Analysis  of  Metamodels 

Residual  analysis  is  intended  to  check  whether  the 
residuals  of  metamodels  are  identically  and  independently 
distributed  (Ilb) ,  which  is  an  important  assumption  of 
linear  regression  analysis. 

The  residual  analysis  in  the  study  was  performed  with 
the  SAS  commercial  statistical  software  package.  The  ARIMA 
process  in  SAS  was  employed  to  test  whether  the  correlation 
of  the  residuals  is  within  the  95%  confidence  interval.  The 
results  are  summarized  as  follows: 

Metamodel  Correlation  95%  Confidence  Interval 


‘^a 

0.17 

0.19 

c  ^ 

0.14 

0.22 

C  ^ 

'“D 

0.20 

0.22 

r  2 

0 . 02 

0.15 

The  results  show  that  the  correlations  are  all  within 
the  95%  confidence  intervals,  indicating  that  the  residuals 
of  each  metamodel  are  identically  and  independently 
distributed.  Therefore,  the  metamodels  developed  in  the 
study  do  not  violate  the  IID  assumption  of  linear 
regression  analysis. 
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